?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
A freeze group can now be allowed to move rigidly in some dimension by using "freezed...
[gromacs/rigid-bodies.git] / src / mdlib / genborn_allvsall_sse2_double.c
blobc91d05a06b9720a5bbd4f5314e7e7f5b5ebaa78f
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 ( (defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || defined(GMX_SSE2)) && defined(GMX_DOUBLE) )
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
2513