src/gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_double/nb_kernel_ElecEw_VdwCSTab_GeomP1P1_avx_128_fma_double.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*
  36  * Note: this file was generated by the GROMACS avx_128_fma_double kernel generator.
  37  */
  38 #include "config.h"
  39
  40 #include <math.h>
  41
  42 #include "../nb_kernel.h"
  43 #include "gromacs/legacyheaders/types/simple.h"
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/legacyheaders/nrnb.h"
  46
  47 #include "gromacs/simd/math_x86_avx_128_fma_double.h"
  48 #include "kernelutil_x86_avx_128_fma_double.h"
  49
  50 /*
  51  * Gromacs nonbonded kernel:   nb_kernel_ElecEw_VdwCSTab_GeomP1P1_VF_avx_128_fma_double
  52  * Electrostatics interaction: Ewald
  53  * VdW interaction:            CubicSplineTable
  54  * Geometry:                   Particle-Particle
  55  * Calculate force/pot:        PotentialAndForce
  56  */
  57 void
  58 nb_kernel_ElecEw_VdwCSTab_GeomP1P1_VF_avx_128_fma_double
  59                     (t_nblist                    * gmx_restrict       nlist,
  60                      rvec                        * gmx_restrict          xx,
  61                      rvec                        * gmx_restrict          ff,
  62                      t_forcerec                  * gmx_restrict          fr,
  63                      t_mdatoms                   * gmx_restrict     mdatoms,
  64                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
  65                      t_nrnb                      * gmx_restrict        nrnb)
  66 {
  67     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
  68      * just 0 for non-waters.
  69      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
  70      * jnr indices corresponding to data put in the four positions in the SIMD register.
  71      */
  72     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
  73     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
  74     int              jnrA,jnrB;
  75     int              j_coord_offsetA,j_coord_offsetB;
  76     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
  77     real             rcutoff_scalar;
  78     real             *shiftvec,*fshift,*x,*f;
  79     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
  80     int              vdwioffset0;
  81     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
  82     int              vdwjidx0A,vdwjidx0B;
  83     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
  84     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
  85     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
  86     real             *charge;
  87     int              nvdwtype;
  88     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
  89     int              *vdwtype;
  90     real             *vdwparam;
  91     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
  92     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
  93     __m128i          vfitab;
  94     __m128i          ifour       = _mm_set1_epi32(4);
  95     __m128d          rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF,twovfeps;
  96     real             *vftab;
  97     __m128i          ewitab;
  98     __m128d          ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
  99     real             *ewtab;
 100     __m128d          dummy_mask,cutoff_mask;
 101     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 102     __m128d          one     = _mm_set1_pd(1.0);
 103     __m128d          two     = _mm_set1_pd(2.0);
 104     x                = xx[0];
 105     f                = ff[0];
 106
 107     nri              = nlist->nri;
 108     iinr             = nlist->iinr;
 109     jindex           = nlist->jindex;
 110     jjnr             = nlist->jjnr;
 111     shiftidx         = nlist->shift;
 112     gid              = nlist->gid;
 113     shiftvec         = fr->shift_vec[0];
 114     fshift           = fr->fshift[0];
 115     facel            = _mm_set1_pd(fr->epsfac);
 116     charge           = mdatoms->chargeA;
 117     nvdwtype         = fr->ntype;
 118     vdwparam         = fr->nbfp;
 119     vdwtype          = mdatoms->typeA;
 120
 121     vftab            = kernel_data->table_vdw->data;
 122     vftabscale       = _mm_set1_pd(kernel_data->table_vdw->scale);
 123
 124     sh_ewald         = _mm_set1_pd(fr->ic->sh_ewald);
 125     ewtab            = fr->ic->tabq_coul_FDV0;
 126     ewtabscale       = _mm_set1_pd(fr->ic->tabq_scale);
 127     ewtabhalfspace   = _mm_set1_pd(0.5/fr->ic->tabq_scale);
 128
 129     /* Avoid stupid compiler warnings */
 130     jnrA = jnrB = 0;
 131     j_coord_offsetA = 0;
 132     j_coord_offsetB = 0;
 133
 134     outeriter        = 0;
 135     inneriter        = 0;
 136
 137     /* Start outer loop over neighborlists */
 138     for(iidx=0; iidx<nri; iidx++)
 139     {
 140         /* Load shift vector for this list */
 141         i_shift_offset   = DIM*shiftidx[iidx];
 142
 143         /* Load limits for loop over neighbors */
 144         j_index_start    = jindex[iidx];
 145         j_index_end      = jindex[iidx+1];
 146
 147         /* Get outer coordinate index */
 148         inr              = iinr[iidx];
 149         i_coord_offset   = DIM*inr;
 150
 151         /* Load i particle coords and add shift vector */
 152         gmx_mm_load_shift_and_1rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,&ix0,&iy0,&iz0);
 153
 154         fix0             = _mm_setzero_pd();
 155         fiy0             = _mm_setzero_pd();
 156         fiz0             = _mm_setzero_pd();
 157
 158         /* Load parameters for i particles */
 159         iq0              = _mm_mul_pd(facel,_mm_load1_pd(charge+inr+0));
 160         vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 161
 162         /* Reset potential sums */
 163         velecsum         = _mm_setzero_pd();
 164         vvdwsum          = _mm_setzero_pd();
 165
 166         /* Start inner kernel loop */
 167         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
 168         {
 169
 170             /* Get j neighbor index, and coordinate index */
 171             jnrA             = jjnr[jidx];
 172             jnrB             = jjnr[jidx+1];
 173             j_coord_offsetA  = DIM*jnrA;
 174             j_coord_offsetB  = DIM*jnrB;
 175
 176             /* load j atom coordinates */
 177             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 178                                               &jx0,&jy0,&jz0);
 179
 180             /* Calculate displacement vector */
 181             dx00             = _mm_sub_pd(ix0,jx0);
 182             dy00             = _mm_sub_pd(iy0,jy0);
 183             dz00             = _mm_sub_pd(iz0,jz0);
 184
 185             /* Calculate squared distance and things based on it */
 186             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 187
 188             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 189
 190             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 191
 192             /* Load parameters for j particles */
 193             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
 194             vdwjidx0A        = 2*vdwtype[jnrA+0];
 195             vdwjidx0B        = 2*vdwtype[jnrB+0];
 196
 197             /**************************
 198              * CALCULATE INTERACTIONS *
 199              **************************/
 200
 201             r00              = _mm_mul_pd(rsq00,rinv00);
 202
 203             /* Compute parameters for interactions between i and j atoms */
 204             qq00             = _mm_mul_pd(iq0,jq0);
 205             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
 206                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
 207
 208             /* Calculate table index by multiplying r with table scale and truncate to integer */
 209             rt               = _mm_mul_pd(r00,vftabscale);
 210             vfitab           = _mm_cvttpd_epi32(rt);
 211 #ifdef __XOP__
 212             vfeps            = _mm_frcz_pd(rt);
 213 #else
 214             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
 215 #endif
 216             twovfeps         = _mm_add_pd(vfeps,vfeps);
 217             vfitab           = _mm_slli_epi32(vfitab,3);
 218
 219             /* EWALD ELECTROSTATICS */
 220
 221             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 222             ewrt             = _mm_mul_pd(r00,ewtabscale);
 223             ewitab           = _mm_cvttpd_epi32(ewrt);
 224 #ifdef __XOP__
 225             eweps            = _mm_frcz_pd(ewrt);
 226 #else
 227             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 228 #endif
 229             twoeweps         = _mm_add_pd(eweps,eweps);
 230             ewitab           = _mm_slli_epi32(ewitab,2);
 231             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 232             ewtabD           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
 233             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 234             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 235             ewtabFn          = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
 236             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 237             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 238             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 239             velec            = _mm_mul_pd(qq00,_mm_sub_pd(rinv00,velec));
 240             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
 241
 242             /* CUBIC SPLINE TABLE DISPERSION */
 243             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
 244             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
 245             GMX_MM_TRANSPOSE2_PD(Y,F);
 246             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
 247             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
 248             GMX_MM_TRANSPOSE2_PD(G,H);
 249             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
 250             VV               = _mm_macc_pd(vfeps,Fp,Y);
 251             vvdw6            = _mm_mul_pd(c6_00,VV);
 252             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
 253             fvdw6            = _mm_mul_pd(c6_00,FF);
 254
 255             /* CUBIC SPLINE TABLE REPULSION */
 256             vfitab           = _mm_add_epi32(vfitab,ifour);
 257             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
 258             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
 259             GMX_MM_TRANSPOSE2_PD(Y,F);
 260             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
 261             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
 262             GMX_MM_TRANSPOSE2_PD(G,H);
 263             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
 264             VV               = _mm_macc_pd(vfeps,Fp,Y);
 265             vvdw12           = _mm_mul_pd(c12_00,VV);
 266             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
 267             fvdw12           = _mm_mul_pd(c12_00,FF);
 268             vvdw             = _mm_add_pd(vvdw12,vvdw6);
 269             fvdw             = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
 270
 271             /* Update potential sum for this i atom from the interaction with this j atom. */
 272             velecsum         = _mm_add_pd(velecsum,velec);
 273             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 274
 275             fscal            = _mm_add_pd(felec,fvdw);
 276
 277             /* Update vectorial force */
 278             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 279             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 280             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 281
 282             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,
 283                                                    _mm_mul_pd(dx00,fscal),
 284                                                    _mm_mul_pd(dy00,fscal),
 285                                                    _mm_mul_pd(dz00,fscal));
 286
 287             /* Inner loop uses 78 flops */
 288         }
 289
 290         if(jidx<j_index_end)
 291         {
 292
 293             jnrA             = jjnr[jidx];
 294             j_coord_offsetA  = DIM*jnrA;
 295
 296             /* load j atom coordinates */
 297             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
 298                                               &jx0,&jy0,&jz0);
 299
 300             /* Calculate displacement vector */
 301             dx00             = _mm_sub_pd(ix0,jx0);
 302             dy00             = _mm_sub_pd(iy0,jy0);
 303             dz00             = _mm_sub_pd(iz0,jz0);
 304
 305             /* Calculate squared distance and things based on it */
 306             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 307
 308             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 309
 310             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 311
 312             /* Load parameters for j particles */
 313             jq0              = _mm_load_sd(charge+jnrA+0);
 314             vdwjidx0A        = 2*vdwtype[jnrA+0];
 315
 316             /**************************
 317              * CALCULATE INTERACTIONS *
 318              **************************/
 319
 320             r00              = _mm_mul_pd(rsq00,rinv00);
 321
 322             /* Compute parameters for interactions between i and j atoms */
 323             qq00             = _mm_mul_pd(iq0,jq0);
 324             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
 325
 326             /* Calculate table index by multiplying r with table scale and truncate to integer */
 327             rt               = _mm_mul_pd(r00,vftabscale);
 328             vfitab           = _mm_cvttpd_epi32(rt);
 329 #ifdef __XOP__
 330             vfeps            = _mm_frcz_pd(rt);
 331 #else
 332             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
 333 #endif
 334             twovfeps         = _mm_add_pd(vfeps,vfeps);
 335             vfitab           = _mm_slli_epi32(vfitab,3);
 336
 337             /* EWALD ELECTROSTATICS */
 338
 339             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 340             ewrt             = _mm_mul_pd(r00,ewtabscale);
 341             ewitab           = _mm_cvttpd_epi32(ewrt);
 342 #ifdef __XOP__
 343             eweps            = _mm_frcz_pd(ewrt);
 344 #else
 345             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 346 #endif
 347             twoeweps         = _mm_add_pd(eweps,eweps);
 348             ewitab           = _mm_slli_epi32(ewitab,2);
 349             ewtabF           = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
 350             ewtabD           = _mm_setzero_pd();
 351             GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
 352             ewtabV           = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
 353             ewtabFn          = _mm_setzero_pd();
 354             GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
 355             felec            = _mm_macc_pd(eweps,ewtabD,ewtabF);
 356             velec            = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
 357             velec            = _mm_mul_pd(qq00,_mm_sub_pd(rinv00,velec));
 358             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
 359
 360             /* CUBIC SPLINE TABLE DISPERSION */
 361             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
 362             F                = _mm_setzero_pd();
 363             GMX_MM_TRANSPOSE2_PD(Y,F);
 364             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
 365             H                = _mm_setzero_pd();
 366             GMX_MM_TRANSPOSE2_PD(G,H);
 367             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
 368             VV               = _mm_macc_pd(vfeps,Fp,Y);
 369             vvdw6            = _mm_mul_pd(c6_00,VV);
 370             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
 371             fvdw6            = _mm_mul_pd(c6_00,FF);
 372
 373             /* CUBIC SPLINE TABLE REPULSION */
 374             vfitab           = _mm_add_epi32(vfitab,ifour);
 375             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
 376             F                = _mm_setzero_pd();
 377             GMX_MM_TRANSPOSE2_PD(Y,F);
 378             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
 379             H                = _mm_setzero_pd();
 380             GMX_MM_TRANSPOSE2_PD(G,H);
 381             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
 382             VV               = _mm_macc_pd(vfeps,Fp,Y);
 383             vvdw12           = _mm_mul_pd(c12_00,VV);
 384             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
 385             fvdw12           = _mm_mul_pd(c12_00,FF);
 386             vvdw             = _mm_add_pd(vvdw12,vvdw6);
 387             fvdw             = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
 388
 389             /* Update potential sum for this i atom from the interaction with this j atom. */
 390             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
 391             velecsum         = _mm_add_pd(velecsum,velec);
 392             vvdw             = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
 393             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
 394
 395             fscal            = _mm_add_pd(felec,fvdw);
 396
 397             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 398
 399             /* Update vectorial force */
 400             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 401             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 402             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 403
 404             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,
 405                                                    _mm_mul_pd(dx00,fscal),
 406                                                    _mm_mul_pd(dy00,fscal),
 407                                                    _mm_mul_pd(dz00,fscal));
 408
 409             /* Inner loop uses 78 flops */
 410         }
 411
 412         /* End of innermost loop */
 413
 414         gmx_mm_update_iforce_1atom_swizzle_pd(fix0,fiy0,fiz0,
 415                                               f+i_coord_offset,fshift+i_shift_offset);
 416
 417         ggid                        = gid[iidx];
 418         /* Update potential energies */
 419         gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
 420         gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
 421
 422         /* Increment number of inner iterations */
 423         inneriter                  += j_index_end - j_index_start;
 424
 425         /* Outer loop uses 9 flops */
 426     }
 427
 428     /* Increment number of outer iterations */
 429     outeriter        += nri;
 430
 431     /* Update outer/inner flops */
 432
 433     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_VF,outeriter*9 + inneriter*78);
 434 }
 435 /*
 436  * Gromacs nonbonded kernel:   nb_kernel_ElecEw_VdwCSTab_GeomP1P1_F_avx_128_fma_double
 437  * Electrostatics interaction: Ewald
 438  * VdW interaction:            CubicSplineTable
 439  * Geometry:                   Particle-Particle
 440  * Calculate force/pot:        Force
 441  */
 442 void
 443 nb_kernel_ElecEw_VdwCSTab_GeomP1P1_F_avx_128_fma_double
 444                     (t_nblist                    * gmx_restrict       nlist,
 445                      rvec                        * gmx_restrict          xx,
 446                      rvec                        * gmx_restrict          ff,
 447                      t_forcerec                  * gmx_restrict          fr,
 448                      t_mdatoms                   * gmx_restrict     mdatoms,
 449                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
 450                      t_nrnb                      * gmx_restrict        nrnb)
 451 {
 452     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
 453      * just 0 for non-waters.
 454      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
 455      * jnr indices corresponding to data put in the four positions in the SIMD register.
 456      */
 457     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
 458     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
 459     int              jnrA,jnrB;
 460     int              j_coord_offsetA,j_coord_offsetB;
 461     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
 462     real             rcutoff_scalar;
 463     real             *shiftvec,*fshift,*x,*f;
 464     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
 465     int              vdwioffset0;
 466     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
 467     int              vdwjidx0A,vdwjidx0B;
 468     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
 469     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
 470     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
 471     real             *charge;
 472     int              nvdwtype;
 473     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
 474     int              *vdwtype;
 475     real             *vdwparam;
 476     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
 477     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
 478     __m128i          vfitab;
 479     __m128i          ifour       = _mm_set1_epi32(4);
 480     __m128d          rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF,twovfeps;
 481     real             *vftab;
 482     __m128i          ewitab;
 483     __m128d          ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
 484     real             *ewtab;
 485     __m128d          dummy_mask,cutoff_mask;
 486     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 487     __m128d          one     = _mm_set1_pd(1.0);
 488     __m128d          two     = _mm_set1_pd(2.0);
 489     x                = xx[0];
 490     f                = ff[0];
 491
 492     nri              = nlist->nri;
 493     iinr             = nlist->iinr;
 494     jindex           = nlist->jindex;
 495     jjnr             = nlist->jjnr;
 496     shiftidx         = nlist->shift;
 497     gid              = nlist->gid;
 498     shiftvec         = fr->shift_vec[0];
 499     fshift           = fr->fshift[0];
 500     facel            = _mm_set1_pd(fr->epsfac);
 501     charge           = mdatoms->chargeA;
 502     nvdwtype         = fr->ntype;
 503     vdwparam         = fr->nbfp;
 504     vdwtype          = mdatoms->typeA;
 505
 506     vftab            = kernel_data->table_vdw->data;
 507     vftabscale       = _mm_set1_pd(kernel_data->table_vdw->scale);
 508
 509     sh_ewald         = _mm_set1_pd(fr->ic->sh_ewald);
 510     ewtab            = fr->ic->tabq_coul_F;
 511     ewtabscale       = _mm_set1_pd(fr->ic->tabq_scale);
 512     ewtabhalfspace   = _mm_set1_pd(0.5/fr->ic->tabq_scale);
 513
 514     /* Avoid stupid compiler warnings */
 515     jnrA = jnrB = 0;
 516     j_coord_offsetA = 0;
 517     j_coord_offsetB = 0;
 518
 519     outeriter        = 0;
 520     inneriter        = 0;
 521
 522     /* Start outer loop over neighborlists */
 523     for(iidx=0; iidx<nri; iidx++)
 524     {
 525         /* Load shift vector for this list */
 526         i_shift_offset   = DIM*shiftidx[iidx];
 527
 528         /* Load limits for loop over neighbors */
 529         j_index_start    = jindex[iidx];
 530         j_index_end      = jindex[iidx+1];
 531
 532         /* Get outer coordinate index */
 533         inr              = iinr[iidx];
 534         i_coord_offset   = DIM*inr;
 535
 536         /* Load i particle coords and add shift vector */
 537         gmx_mm_load_shift_and_1rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,&ix0,&iy0,&iz0);
 538
 539         fix0             = _mm_setzero_pd();
 540         fiy0             = _mm_setzero_pd();
 541         fiz0             = _mm_setzero_pd();
 542
 543         /* Load parameters for i particles */
 544         iq0              = _mm_mul_pd(facel,_mm_load1_pd(charge+inr+0));
 545         vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
 546
 547         /* Start inner kernel loop */
 548         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
 549         {
 550
 551             /* Get j neighbor index, and coordinate index */
 552             jnrA             = jjnr[jidx];
 553             jnrB             = jjnr[jidx+1];
 554             j_coord_offsetA  = DIM*jnrA;
 555             j_coord_offsetB  = DIM*jnrB;
 556
 557             /* load j atom coordinates */
 558             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
 559                                               &jx0,&jy0,&jz0);
 560
 561             /* Calculate displacement vector */
 562             dx00             = _mm_sub_pd(ix0,jx0);
 563             dy00             = _mm_sub_pd(iy0,jy0);
 564             dz00             = _mm_sub_pd(iz0,jz0);
 565
 566             /* Calculate squared distance and things based on it */
 567             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 568
 569             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 570
 571             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 572
 573             /* Load parameters for j particles */
 574             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
 575             vdwjidx0A        = 2*vdwtype[jnrA+0];
 576             vdwjidx0B        = 2*vdwtype[jnrB+0];
 577
 578             /**************************
 579              * CALCULATE INTERACTIONS *
 580              **************************/
 581
 582             r00              = _mm_mul_pd(rsq00,rinv00);
 583
 584             /* Compute parameters for interactions between i and j atoms */
 585             qq00             = _mm_mul_pd(iq0,jq0);
 586             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
 587                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
 588
 589             /* Calculate table index by multiplying r with table scale and truncate to integer */
 590             rt               = _mm_mul_pd(r00,vftabscale);
 591             vfitab           = _mm_cvttpd_epi32(rt);
 592 #ifdef __XOP__
 593             vfeps            = _mm_frcz_pd(rt);
 594 #else
 595             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
 596 #endif
 597             twovfeps         = _mm_add_pd(vfeps,vfeps);
 598             vfitab           = _mm_slli_epi32(vfitab,3);
 599
 600             /* EWALD ELECTROSTATICS */
 601
 602             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 603             ewrt             = _mm_mul_pd(r00,ewtabscale);
 604             ewitab           = _mm_cvttpd_epi32(ewrt);
 605 #ifdef __XOP__
 606             eweps            = _mm_frcz_pd(ewrt);
 607 #else
 608             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 609 #endif
 610             twoeweps         = _mm_add_pd(eweps,eweps);
 611             gmx_mm_load_2pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),ewtab+_mm_extract_epi32(ewitab,1),
 612                                          &ewtabF,&ewtabFn);
 613             felec            = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
 614             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
 615
 616             /* CUBIC SPLINE TABLE DISPERSION */
 617             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
 618             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
 619             GMX_MM_TRANSPOSE2_PD(Y,F);
 620             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
 621             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
 622             GMX_MM_TRANSPOSE2_PD(G,H);
 623             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
 624             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
 625             fvdw6            = _mm_mul_pd(c6_00,FF);
 626
 627             /* CUBIC SPLINE TABLE REPULSION */
 628             vfitab           = _mm_add_epi32(vfitab,ifour);
 629             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
 630             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
 631             GMX_MM_TRANSPOSE2_PD(Y,F);
 632             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
 633             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
 634             GMX_MM_TRANSPOSE2_PD(G,H);
 635             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
 636             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
 637             fvdw12           = _mm_mul_pd(c12_00,FF);
 638             fvdw             = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
 639
 640             fscal            = _mm_add_pd(felec,fvdw);
 641
 642             /* Update vectorial force */
 643             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 644             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 645             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 646
 647             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,
 648                                                    _mm_mul_pd(dx00,fscal),
 649                                                    _mm_mul_pd(dy00,fscal),
 650                                                    _mm_mul_pd(dz00,fscal));
 651
 652             /* Inner loop uses 65 flops */
 653         }
 654
 655         if(jidx<j_index_end)
 656         {
 657
 658             jnrA             = jjnr[jidx];
 659             j_coord_offsetA  = DIM*jnrA;
 660
 661             /* load j atom coordinates */
 662             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
 663                                               &jx0,&jy0,&jz0);
 664
 665             /* Calculate displacement vector */
 666             dx00             = _mm_sub_pd(ix0,jx0);
 667             dy00             = _mm_sub_pd(iy0,jy0);
 668             dz00             = _mm_sub_pd(iz0,jz0);
 669
 670             /* Calculate squared distance and things based on it */
 671             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
 672
 673             rinv00           = gmx_mm_invsqrt_pd(rsq00);
 674
 675             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
 676
 677             /* Load parameters for j particles */
 678             jq0              = _mm_load_sd(charge+jnrA+0);
 679             vdwjidx0A        = 2*vdwtype[jnrA+0];
 680
 681             /**************************
 682              * CALCULATE INTERACTIONS *
 683              **************************/
 684
 685             r00              = _mm_mul_pd(rsq00,rinv00);
 686
 687             /* Compute parameters for interactions between i and j atoms */
 688             qq00             = _mm_mul_pd(iq0,jq0);
 689             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
 690
 691             /* Calculate table index by multiplying r with table scale and truncate to integer */
 692             rt               = _mm_mul_pd(r00,vftabscale);
 693             vfitab           = _mm_cvttpd_epi32(rt);
 694 #ifdef __XOP__
 695             vfeps            = _mm_frcz_pd(rt);
 696 #else
 697             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
 698 #endif
 699             twovfeps         = _mm_add_pd(vfeps,vfeps);
 700             vfitab           = _mm_slli_epi32(vfitab,3);
 701
 702             /* EWALD ELECTROSTATICS */
 703
 704             /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
 705             ewrt             = _mm_mul_pd(r00,ewtabscale);
 706             ewitab           = _mm_cvttpd_epi32(ewrt);
 707 #ifdef __XOP__
 708             eweps            = _mm_frcz_pd(ewrt);
 709 #else
 710             eweps            = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
 711 #endif
 712             twoeweps         = _mm_add_pd(eweps,eweps);
 713             gmx_mm_load_1pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
 714             felec            = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
 715             felec            = _mm_mul_pd(_mm_mul_pd(qq00,rinv00),_mm_sub_pd(rinvsq00,felec));
 716
 717             /* CUBIC SPLINE TABLE DISPERSION */
 718             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
 719             F                = _mm_setzero_pd();
 720             GMX_MM_TRANSPOSE2_PD(Y,F);
 721             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
 722             H                = _mm_setzero_pd();
 723             GMX_MM_TRANSPOSE2_PD(G,H);
 724             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
 725             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
 726             fvdw6            = _mm_mul_pd(c6_00,FF);
 727
 728             /* CUBIC SPLINE TABLE REPULSION */
 729             vfitab           = _mm_add_epi32(vfitab,ifour);
 730             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
 731             F                = _mm_setzero_pd();
 732             GMX_MM_TRANSPOSE2_PD(Y,F);
 733             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
 734             H                = _mm_setzero_pd();
 735             GMX_MM_TRANSPOSE2_PD(G,H);
 736             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
 737             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
 738             fvdw12           = _mm_mul_pd(c12_00,FF);
 739             fvdw             = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
 740
 741             fscal            = _mm_add_pd(felec,fvdw);
 742
 743             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
 744
 745             /* Update vectorial force */
 746             fix0             = _mm_macc_pd(dx00,fscal,fix0);
 747             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
 748             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
 749
 750             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,
 751                                                    _mm_mul_pd(dx00,fscal),
 752                                                    _mm_mul_pd(dy00,fscal),
 753                                                    _mm_mul_pd(dz00,fscal));
 754
 755             /* Inner loop uses 65 flops */
 756         }
 757
 758         /* End of innermost loop */
 759
 760         gmx_mm_update_iforce_1atom_swizzle_pd(fix0,fiy0,fiz0,
 761                                               f+i_coord_offset,fshift+i_shift_offset);
 762
 763         /* Increment number of inner iterations */
 764         inneriter                  += j_index_end - j_index_start;
 765
 766         /* Outer loop uses 7 flops */
 767     }
 768
 769     /* Increment number of outer iterations */
 770     outeriter        += nri;
 771
 772     /* Update outer/inner flops */
 773
 774     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_F,outeriter*7 + inneriter*65);
 775 }