1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2008, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Groningen Machine for Chemical Simulation
53 #include "gmx_fatal.h"
54 #include "mtop_util.h"
57 #include "gromacs/utility/gmxmpi.h"
59 /* Only compile this file if SSE2 intrinsics are available */
60 #if 0 && defined (GMX_X86_SSE2)
61 #include <gmx_sse2_double.h>
62 #include <emmintrin.h>
64 #include "genborn_sse2_double.h"
67 calc_gb_rad_still_sse2_double(t_commrec *cr, t_forcerec *fr,
68 int natoms, gmx_localtop_t *top,
69 double *x, t_nblist *nl,
72 int i, k, n, ii, is3, ii3, nj0, nj1, offset;
73 int jnrA, jnrB, j3A, j3B;
90 __m128d rsq, rinv, rinv2, rinv4, rinv6;
91 __m128d ratio, gpi, rai, raj, vai, vaj, rvdw;
92 __m128d ccf, dccf, theta, cosq, term, sinq, res, prod, prod_ai, tmp;
93 __m128d mask, icf4, icf6, mask_cmp;
95 const __m128d half = _mm_set1_pd(0.5);
96 const __m128d three = _mm_set1_pd(3.0);
97 const __m128d one = _mm_set1_pd(1.0);
98 const __m128d two = _mm_set1_pd(2.0);
99 const __m128d zero = _mm_set1_pd(0.0);
100 const __m128d four = _mm_set1_pd(4.0);
102 const __m128d still_p5inv = _mm_set1_pd(STILL_P5INV);
103 const __m128d still_pip5 = _mm_set1_pd(STILL_PIP5);
104 const __m128d still_p4 = _mm_set1_pd(STILL_P4);
106 factor = 0.5 * ONE_4PI_EPS0;
108 gb_radius = born->gb_radius;
110 work = born->gpol_still_work;
112 shiftvec = fr->shift_vec[0];
116 jx = _mm_setzero_pd();
117 jy = _mm_setzero_pd();
118 jz = _mm_setzero_pd();
122 for (i = 0; i < natoms; i++)
127 for (i = 0; i < nl->nri; i++)
131 is3 = 3*nl->shift[i];
133 shY = shiftvec[is3+1];
134 shZ = shiftvec[is3+2];
136 nj1 = nl->jindex[i+1];
138 ix = _mm_set1_pd(shX+x[ii3+0]);
139 iy = _mm_set1_pd(shY+x[ii3+1]);
140 iz = _mm_set1_pd(shZ+x[ii3+2]);
143 /* Polarization energy for atom ai */
144 gpi = _mm_setzero_pd();
146 rai = _mm_load1_pd(gb_radius+ii);
147 prod_ai = _mm_set1_pd(STILL_P4*vsolv[ii]);
149 for (k = nj0; k < nj1-1; k += 2)
157 GMX_MM_LOAD_1RVEC_2POINTERS_PD(x+j3A, x+j3B, jx, jy, jz);
159 GMX_MM_LOAD_2VALUES_PD(gb_radius+jnrA, gb_radius+jnrB, raj);
160 GMX_MM_LOAD_2VALUES_PD(vsolv+jnrA, vsolv+jnrB, vaj);
162 dx = _mm_sub_pd(ix, jx);
163 dy = _mm_sub_pd(iy, jy);
164 dz = _mm_sub_pd(iz, jz);
166 rsq = gmx_mm_calc_rsq_pd(dx, dy, dz);
167 rinv = gmx_mm_invsqrt_pd(rsq);
168 rinv2 = _mm_mul_pd(rinv, rinv);
169 rinv4 = _mm_mul_pd(rinv2, rinv2);
170 rinv6 = _mm_mul_pd(rinv4, rinv2);
172 rvdw = _mm_add_pd(rai, raj);
173 ratio = _mm_mul_pd(rsq, gmx_mm_inv_pd( _mm_mul_pd(rvdw, rvdw)));
175 mask_cmp = _mm_cmple_pd(ratio, still_p5inv);
177 /* gmx_mm_sincos_pd() is quite expensive, so avoid calculating it if we can! */
178 if (0 == _mm_movemask_pd(mask_cmp) )
180 /* if ratio>still_p5inv for ALL elements */
182 dccf = _mm_setzero_pd();
186 ratio = _mm_min_pd(ratio, still_p5inv);
187 theta = _mm_mul_pd(ratio, still_pip5);
188 gmx_mm_sincos_pd(theta, &sinq, &cosq);
189 term = _mm_mul_pd(half, _mm_sub_pd(one, cosq));
190 ccf = _mm_mul_pd(term, term);
191 dccf = _mm_mul_pd(_mm_mul_pd(two, term),
192 _mm_mul_pd(sinq, theta));
195 prod = _mm_mul_pd(still_p4, vaj);
196 icf4 = _mm_mul_pd(ccf, rinv4);
197 icf6 = _mm_mul_pd( _mm_sub_pd( _mm_mul_pd(four, ccf), dccf), rinv6);
199 GMX_MM_INCREMENT_2VALUES_PD(work+jnrA, work+jnrB, _mm_mul_pd(prod_ai, icf4));
201 gpi = _mm_add_pd(gpi, _mm_mul_pd(prod, icf4) );
203 _mm_store_pd(dadx, _mm_mul_pd(prod, icf6));
205 _mm_store_pd(dadx, _mm_mul_pd(prod_ai, icf6));
215 GMX_MM_LOAD_1RVEC_1POINTER_PD(x+j3A, jx, jy, jz);
217 GMX_MM_LOAD_1VALUE_PD(gb_radius+jnrA, raj);
218 GMX_MM_LOAD_1VALUE_PD(vsolv+jnrA, vaj);
220 dx = _mm_sub_sd(ix, jx);
221 dy = _mm_sub_sd(iy, jy);
222 dz = _mm_sub_sd(iz, jz);
224 rsq = gmx_mm_calc_rsq_pd(dx, dy, dz);
225 rinv = gmx_mm_invsqrt_pd(rsq);
226 rinv2 = _mm_mul_sd(rinv, rinv);
227 rinv4 = _mm_mul_sd(rinv2, rinv2);
228 rinv6 = _mm_mul_sd(rinv4, rinv2);
230 rvdw = _mm_add_sd(rai, raj);
231 ratio = _mm_mul_sd(rsq, gmx_mm_inv_pd( _mm_mul_pd(rvdw, rvdw)));
233 mask_cmp = _mm_cmple_sd(ratio, still_p5inv);
235 /* gmx_mm_sincos_pd() is quite expensive, so avoid calculating it if we can! */
236 if (0 == _mm_movemask_pd(mask_cmp) )
238 /* if ratio>still_p5inv for ALL elements */
240 dccf = _mm_setzero_pd();
244 ratio = _mm_min_sd(ratio, still_p5inv);
245 theta = _mm_mul_sd(ratio, still_pip5);
246 gmx_mm_sincos_pd(theta, &sinq, &cosq);
247 term = _mm_mul_sd(half, _mm_sub_sd(one, cosq));
248 ccf = _mm_mul_sd(term, term);
249 dccf = _mm_mul_sd(_mm_mul_sd(two, term),
250 _mm_mul_sd(sinq, theta));
253 prod = _mm_mul_sd(still_p4, vaj);
254 icf4 = _mm_mul_sd(ccf, rinv4);
255 icf6 = _mm_mul_sd( _mm_sub_sd( _mm_mul_sd(four, ccf), dccf), rinv6);
257 GMX_MM_INCREMENT_1VALUE_PD(work+jnrA, _mm_mul_sd(prod_ai, icf4));
259 gpi = _mm_add_sd(gpi, _mm_mul_sd(prod, icf4) );
261 _mm_store_pd(dadx, _mm_mul_pd(prod, icf6));
263 _mm_store_pd(dadx, _mm_mul_pd(prod_ai, icf6));
266 gmx_mm_update_1pot_pd(gpi, work+ii);
269 /* Sum up the polarization energy from other nodes */
272 gmx_sum(natoms, work, cr);
274 else if (DOMAINDECOMP(cr))
276 dd_atom_sum_real(cr->dd, work);
279 /* Compute the radii */
280 for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
282 if (born->use[i] != 0)
284 gpi_ai = born->gpol[i] + work[i]; /* add gpi to the initial pol energy gpi_ai*/
285 gpi2 = gpi_ai * gpi_ai;
286 born->bRad[i] = factor*gmx_invsqrt(gpi2);
287 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
291 /* Extra (local) communication required for DD */
292 if (DOMAINDECOMP(cr))
294 dd_atom_spread_real(cr->dd, born->bRad);
295 dd_atom_spread_real(cr->dd, fr->invsqrta);
303 calc_gb_rad_hct_obc_sse2_double(t_commrec *cr, t_forcerec * fr, int natoms, gmx_localtop_t *top,
304 double *x, t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md, int gb_algorithm)
306 int i, ai, k, n, ii, ii3, is3, nj0, nj1, at0, at1, offset;
309 double shX, shY, shZ;
310 double rr, rr_inv, rr_inv2, sum_tmp, sum, sum2, sum3, gbr;
311 double sum_ai2, sum_ai3, tsum, tchain, doffset;
320 __m128d ix, iy, iz, jx, jy, jz;
321 __m128d dx, dy, dz, t1, t2, t3, t4;
322 __m128d rsq, rinv, r;
323 __m128d rai, rai_inv, raj, raj_inv, rai_inv2, sk, sk2, lij, dlij, duij;
324 __m128d uij, lij2, uij2, lij3, uij3, diff2;
325 __m128d lij_inv, sk2_inv, prod, log_term, tmp, tmp_sum;
326 __m128d sum_ai, tmp_ai, sk_ai, sk_aj, sk2_ai, sk2_aj, sk2_rinv;
327 __m128d dadx1, dadx2;
330 __m128d obc_mask1, obc_mask2, obc_mask3;
332 __m128d oneeighth = _mm_set1_pd(0.125);
333 __m128d onefourth = _mm_set1_pd(0.25);
335 const __m128d half = _mm_set1_pd(0.5);
336 const __m128d three = _mm_set1_pd(3.0);
337 const __m128d one = _mm_set1_pd(1.0);
338 const __m128d two = _mm_set1_pd(2.0);
339 const __m128d zero = _mm_set1_pd(0.0);
340 const __m128d neg = _mm_set1_pd(-1.0);
342 /* Set the dielectric offset */
343 doffset = born->gb_doffset;
344 gb_radius = born->gb_radius;
345 obc_param = born->param;
346 work = born->gpol_hct_work;
349 shiftvec = fr->shift_vec[0];
351 jx = _mm_setzero_pd();
352 jy = _mm_setzero_pd();
353 jz = _mm_setzero_pd();
357 for (i = 0; i < born->nr; i++)
362 for (i = 0; i < nl->nri; i++)
366 is3 = 3*nl->shift[i];
368 shY = shiftvec[is3+1];
369 shZ = shiftvec[is3+2];
371 nj1 = nl->jindex[i+1];
373 ix = _mm_set1_pd(shX+x[ii3+0]);
374 iy = _mm_set1_pd(shY+x[ii3+1]);
375 iz = _mm_set1_pd(shZ+x[ii3+2]);
377 rai = _mm_load1_pd(gb_radius+ii);
378 rai_inv = gmx_mm_inv_pd(rai);
380 sum_ai = _mm_setzero_pd();
382 sk_ai = _mm_load1_pd(born->param+ii);
383 sk2_ai = _mm_mul_pd(sk_ai, sk_ai);
385 for (k = nj0; k < nj1-1; k += 2)
393 GMX_MM_LOAD_1RVEC_2POINTERS_PD(x+j3A, x+j3B, jx, jy, jz);
394 GMX_MM_LOAD_2VALUES_PD(gb_radius+jnrA, gb_radius+jnrB, raj);
395 GMX_MM_LOAD_2VALUES_PD(obc_param+jnrA, obc_param+jnrB, sk_aj);
397 dx = _mm_sub_pd(ix, jx);
398 dy = _mm_sub_pd(iy, jy);
399 dz = _mm_sub_pd(iz, jz);
401 rsq = gmx_mm_calc_rsq_pd(dx, dy, dz);
403 rinv = gmx_mm_invsqrt_pd(rsq);
404 r = _mm_mul_pd(rsq, rinv);
406 /* Compute raj_inv aj1-4 */
407 raj_inv = gmx_mm_inv_pd(raj);
409 /* Evaluate influence of atom aj -> ai */
410 t1 = _mm_add_pd(r, sk_aj);
411 t2 = _mm_sub_pd(r, sk_aj);
412 t3 = _mm_sub_pd(sk_aj, r);
413 obc_mask1 = _mm_cmplt_pd(rai, t1);
414 obc_mask2 = _mm_cmplt_pd(rai, t2);
415 obc_mask3 = _mm_cmplt_pd(rai, t3);
417 uij = gmx_mm_inv_pd(t1);
418 lij = _mm_or_pd( _mm_and_pd(obc_mask2, gmx_mm_inv_pd(t2)),
419 _mm_andnot_pd(obc_mask2, rai_inv));
420 dlij = _mm_and_pd(one, obc_mask2);
421 uij2 = _mm_mul_pd(uij, uij);
422 uij3 = _mm_mul_pd(uij2, uij);
423 lij2 = _mm_mul_pd(lij, lij);
424 lij3 = _mm_mul_pd(lij2, lij);
426 diff2 = _mm_sub_pd(uij2, lij2);
427 lij_inv = gmx_mm_invsqrt_pd(lij2);
428 sk2_aj = _mm_mul_pd(sk_aj, sk_aj);
429 sk2_rinv = _mm_mul_pd(sk2_aj, rinv);
430 prod = _mm_mul_pd(onefourth, sk2_rinv);
432 logterm = gmx_mm_log_pd(_mm_mul_pd(uij, lij_inv));
434 t1 = _mm_sub_pd(lij, uij);
435 t2 = _mm_mul_pd(diff2,
436 _mm_sub_pd(_mm_mul_pd(onefourth, r),
438 t3 = _mm_mul_pd(half, _mm_mul_pd(rinv, logterm));
439 t1 = _mm_add_pd(t1, _mm_add_pd(t2, t3));
440 t4 = _mm_mul_pd(two, _mm_sub_pd(rai_inv, lij));
441 t4 = _mm_and_pd(t4, obc_mask3);
442 t1 = _mm_mul_pd(half, _mm_add_pd(t1, t4));
444 sum_ai = _mm_add_pd(sum_ai, _mm_and_pd(t1, obc_mask1) );
446 t1 = _mm_add_pd(_mm_mul_pd(half, lij2),
447 _mm_mul_pd(prod, lij3));
449 _mm_mul_pd(onefourth,
450 _mm_add_pd(_mm_mul_pd(lij, rinv),
451 _mm_mul_pd(lij3, r))));
452 t2 = _mm_mul_pd(onefourth,
453 _mm_add_pd(_mm_mul_pd(uij, rinv),
454 _mm_mul_pd(uij3, r)));
456 _mm_add_pd(_mm_mul_pd(half, uij2),
457 _mm_mul_pd(prod, uij3)));
458 t3 = _mm_mul_pd(_mm_mul_pd(onefourth, logterm),
459 _mm_mul_pd(rinv, rinv));
461 _mm_mul_pd(_mm_mul_pd(diff2, oneeighth),
463 _mm_mul_pd(sk2_rinv, rinv))));
464 t1 = _mm_mul_pd(rinv,
465 _mm_add_pd(_mm_mul_pd(dlij, t1),
466 _mm_add_pd(t2, t3)));
468 dadx1 = _mm_and_pd(t1, obc_mask1);
470 /* Evaluate influence of atom ai -> aj */
471 t1 = _mm_add_pd(r, sk_ai);
472 t2 = _mm_sub_pd(r, sk_ai);
473 t3 = _mm_sub_pd(sk_ai, r);
474 obc_mask1 = _mm_cmplt_pd(raj, t1);
475 obc_mask2 = _mm_cmplt_pd(raj, t2);
476 obc_mask3 = _mm_cmplt_pd(raj, t3);
478 uij = gmx_mm_inv_pd(t1);
479 lij = _mm_or_pd( _mm_and_pd(obc_mask2, gmx_mm_inv_pd(t2)),
480 _mm_andnot_pd(obc_mask2, raj_inv));
481 dlij = _mm_and_pd(one, obc_mask2);
482 uij2 = _mm_mul_pd(uij, uij);
483 uij3 = _mm_mul_pd(uij2, uij);
484 lij2 = _mm_mul_pd(lij, lij);
485 lij3 = _mm_mul_pd(lij2, lij);
487 diff2 = _mm_sub_pd(uij2, lij2);
488 lij_inv = gmx_mm_invsqrt_pd(lij2);
489 sk2_rinv = _mm_mul_pd(sk2_ai, rinv);
490 prod = _mm_mul_pd(onefourth, sk2_rinv);
492 logterm = gmx_mm_log_pd(_mm_mul_pd(uij, lij_inv));
494 t1 = _mm_sub_pd(lij, uij);
495 t2 = _mm_mul_pd(diff2,
496 _mm_sub_pd(_mm_mul_pd(onefourth, r),
498 t3 = _mm_mul_pd(half, _mm_mul_pd(rinv, logterm));
499 t1 = _mm_add_pd(t1, _mm_add_pd(t2, t3));
500 t4 = _mm_mul_pd(two, _mm_sub_pd(raj_inv, lij));
501 t4 = _mm_and_pd(t4, obc_mask3);
502 t1 = _mm_mul_pd(half, _mm_add_pd(t1, t4));
504 GMX_MM_INCREMENT_2VALUES_PD(work+jnrA, work+jnrB, _mm_and_pd(t1, obc_mask1));
506 t1 = _mm_add_pd(_mm_mul_pd(half, lij2),
507 _mm_mul_pd(prod, lij3));
509 _mm_mul_pd(onefourth,
510 _mm_add_pd(_mm_mul_pd(lij, rinv),
511 _mm_mul_pd(lij3, r))));
512 t2 = _mm_mul_pd(onefourth,
513 _mm_add_pd(_mm_mul_pd(uij, rinv),
514 _mm_mul_pd(uij3, r)));
516 _mm_add_pd(_mm_mul_pd(half, uij2),
517 _mm_mul_pd(prod, uij3)));
518 t3 = _mm_mul_pd(_mm_mul_pd(onefourth, logterm),
519 _mm_mul_pd(rinv, rinv));
521 _mm_mul_pd(_mm_mul_pd(diff2, oneeighth),
523 _mm_mul_pd(sk2_rinv, rinv))));
524 t1 = _mm_mul_pd(rinv,
525 _mm_add_pd(_mm_mul_pd(dlij, t1),
526 _mm_add_pd(t2, t3)));
528 dadx2 = _mm_and_pd(t1, obc_mask1);
530 _mm_store_pd(dadx, dadx1);
532 _mm_store_pd(dadx, dadx2);
534 } /* end normal inner loop */
542 GMX_MM_LOAD_1RVEC_1POINTER_PD(x+j3A, jx, jy, jz);
543 GMX_MM_LOAD_1VALUE_PD(gb_radius+jnrA, raj);
544 GMX_MM_LOAD_1VALUE_PD(obc_param+jnrA, sk_aj);
546 dx = _mm_sub_sd(ix, jx);
547 dy = _mm_sub_sd(iy, jy);
548 dz = _mm_sub_sd(iz, jz);
550 rsq = gmx_mm_calc_rsq_pd(dx, dy, dz);
552 rinv = gmx_mm_invsqrt_pd(rsq);
553 r = _mm_mul_sd(rsq, rinv);
555 /* Compute raj_inv aj1-4 */
556 raj_inv = gmx_mm_inv_pd(raj);
558 /* Evaluate influence of atom aj -> ai */
559 t1 = _mm_add_sd(r, sk_aj);
560 t2 = _mm_sub_sd(r, sk_aj);
561 t3 = _mm_sub_sd(sk_aj, r);
562 obc_mask1 = _mm_cmplt_sd(rai, t1);
563 obc_mask2 = _mm_cmplt_sd(rai, t2);
564 obc_mask3 = _mm_cmplt_sd(rai, t3);
566 uij = gmx_mm_inv_pd(t1);
567 lij = _mm_or_pd(_mm_and_pd(obc_mask2, gmx_mm_inv_pd(t2)),
568 _mm_andnot_pd(obc_mask2, rai_inv));
569 dlij = _mm_and_pd(one, obc_mask2);
570 uij2 = _mm_mul_sd(uij, uij);
571 uij3 = _mm_mul_sd(uij2, uij);
572 lij2 = _mm_mul_sd(lij, lij);
573 lij3 = _mm_mul_sd(lij2, lij);
575 diff2 = _mm_sub_sd(uij2, lij2);
576 lij_inv = gmx_mm_invsqrt_pd(lij2);
577 sk2_aj = _mm_mul_sd(sk_aj, sk_aj);
578 sk2_rinv = _mm_mul_sd(sk2_aj, rinv);
579 prod = _mm_mul_sd(onefourth, sk2_rinv);
581 logterm = gmx_mm_log_pd(_mm_mul_sd(uij, lij_inv));
583 t1 = _mm_sub_sd(lij, uij);
584 t2 = _mm_mul_sd(diff2,
585 _mm_sub_sd(_mm_mul_pd(onefourth, r),
587 t3 = _mm_mul_sd(half, _mm_mul_sd(rinv, logterm));
588 t1 = _mm_add_sd(t1, _mm_add_sd(t2, t3));
589 t4 = _mm_mul_sd(two, _mm_sub_sd(rai_inv, lij));
590 t4 = _mm_and_pd(t4, obc_mask3);
591 t1 = _mm_mul_sd(half, _mm_add_sd(t1, t4));
593 sum_ai = _mm_add_sd(sum_ai, _mm_and_pd(t1, obc_mask1) );
595 t1 = _mm_add_sd(_mm_mul_sd(half, lij2),
596 _mm_mul_sd(prod, lij3));
598 _mm_mul_sd(onefourth,
599 _mm_add_sd(_mm_mul_sd(lij, rinv),
600 _mm_mul_sd(lij3, r))));
601 t2 = _mm_mul_sd(onefourth,
602 _mm_add_sd(_mm_mul_sd(uij, rinv),
603 _mm_mul_sd(uij3, r)));
605 _mm_add_sd(_mm_mul_sd(half, uij2),
606 _mm_mul_sd(prod, uij3)));
607 t3 = _mm_mul_sd(_mm_mul_sd(onefourth, logterm),
608 _mm_mul_sd(rinv, rinv));
610 _mm_mul_sd(_mm_mul_sd(diff2, oneeighth),
612 _mm_mul_sd(sk2_rinv, rinv))));
613 t1 = _mm_mul_sd(rinv,
614 _mm_add_sd(_mm_mul_sd(dlij, t1),
615 _mm_add_pd(t2, t3)));
617 dadx1 = _mm_and_pd(t1, obc_mask1);
619 /* Evaluate influence of atom ai -> aj */
620 t1 = _mm_add_sd(r, sk_ai);
621 t2 = _mm_sub_sd(r, sk_ai);
622 t3 = _mm_sub_sd(sk_ai, r);
623 obc_mask1 = _mm_cmplt_sd(raj, t1);
624 obc_mask2 = _mm_cmplt_sd(raj, t2);
625 obc_mask3 = _mm_cmplt_sd(raj, t3);
627 uij = gmx_mm_inv_pd(t1);
628 lij = _mm_or_pd( _mm_and_pd(obc_mask2, gmx_mm_inv_pd(t2)),
629 _mm_andnot_pd(obc_mask2, raj_inv));
630 dlij = _mm_and_pd(one, obc_mask2);
631 uij2 = _mm_mul_sd(uij, uij);
632 uij3 = _mm_mul_sd(uij2, uij);
633 lij2 = _mm_mul_sd(lij, lij);
634 lij3 = _mm_mul_sd(lij2, lij);
636 diff2 = _mm_sub_sd(uij2, lij2);
637 lij_inv = gmx_mm_invsqrt_pd(lij2);
638 sk2_rinv = _mm_mul_sd(sk2_ai, rinv);
639 prod = _mm_mul_sd(onefourth, sk2_rinv);
641 logterm = gmx_mm_log_pd(_mm_mul_sd(uij, lij_inv));
643 t1 = _mm_sub_sd(lij, uij);
644 t2 = _mm_mul_sd(diff2,
645 _mm_sub_sd(_mm_mul_sd(onefourth, r),
647 t3 = _mm_mul_sd(half, _mm_mul_sd(rinv, logterm));
648 t1 = _mm_add_sd(t1, _mm_add_sd(t2, t3));
649 t4 = _mm_mul_sd(two, _mm_sub_sd(raj_inv, lij));
650 t4 = _mm_and_pd(t4, obc_mask3);
651 t1 = _mm_mul_sd(half, _mm_add_sd(t1, t4));
653 GMX_MM_INCREMENT_1VALUE_PD(work+jnrA, _mm_and_pd(t1, obc_mask1));
655 t1 = _mm_add_sd(_mm_mul_sd(half, lij2),
656 _mm_mul_sd(prod, lij3));
658 _mm_mul_sd(onefourth,
659 _mm_add_sd(_mm_mul_sd(lij, rinv),
660 _mm_mul_sd(lij3, r))));
661 t2 = _mm_mul_sd(onefourth,
662 _mm_add_sd(_mm_mul_sd(uij, rinv),
663 _mm_mul_sd(uij3, r)));
665 _mm_add_sd(_mm_mul_sd(half, uij2),
666 _mm_mul_sd(prod, uij3)));
667 t3 = _mm_mul_sd(_mm_mul_sd(onefourth, logterm),
668 _mm_mul_sd(rinv, rinv));
670 _mm_mul_sd(_mm_mul_sd(diff2, oneeighth),
672 _mm_mul_sd(sk2_rinv, rinv))));
673 t1 = _mm_mul_sd(rinv,
674 _mm_add_sd(_mm_mul_sd(dlij, t1),
675 _mm_add_sd(t2, t3)));
677 dadx2 = _mm_and_pd(t1, obc_mask1);
679 _mm_store_pd(dadx, dadx1);
681 _mm_store_pd(dadx, dadx2);
684 gmx_mm_update_1pot_pd(sum_ai, work+ii);
688 /* Parallel summations */
691 gmx_sum(natoms, work, cr);
693 else if (DOMAINDECOMP(cr))
695 dd_atom_sum_real(cr->dd, work);
698 if (gb_algorithm == egbHCT)
701 for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
703 if (born->use[i] != 0)
705 rr = top->atomtypes.gb_radius[md->typeA[i]]-doffset;
706 sum = 1.0/rr - work[i];
707 min_rad = rr + doffset;
710 born->bRad[i] = rad > min_rad ? rad : min_rad;
711 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
715 /* Extra communication required for DD */
716 if (DOMAINDECOMP(cr))
718 dd_atom_spread_real(cr->dd, born->bRad);
719 dd_atom_spread_real(cr->dd, fr->invsqrta);
725 for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
727 if (born->use[i] != 0)
729 rr = top->atomtypes.gb_radius[md->typeA[i]];
737 tsum = tanh(born->obc_alpha*sum-born->obc_beta*sum2+born->obc_gamma*sum3);
738 born->bRad[i] = rr_inv - tsum*rr_inv2;
739 born->bRad[i] = 1.0 / born->bRad[i];
741 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
743 tchain = rr * (born->obc_alpha-2*born->obc_beta*sum+3*born->obc_gamma*sum2);
744 born->drobc[i] = (1.0-tsum*tsum)*tchain*rr_inv2;
747 /* Extra (local) communication required for DD */
748 if (DOMAINDECOMP(cr))
750 dd_atom_spread_real(cr->dd, born->bRad);
751 dd_atom_spread_real(cr->dd, fr->invsqrta);
752 dd_atom_spread_real(cr->dd, born->drobc);
763 calc_gb_chainrule_sse2_double(int natoms, t_nblist *nl, double *dadx, double *dvda,
764 double *x, double *f, double *fshift, double *shiftvec,
765 int gb_algorithm, gmx_genborn_t *born, t_mdatoms *md)
767 int i, k, n, ii, jnr, ii3, is3, nj0, nj1, n0, n1;
772 double rbi, shX, shY, shZ;
777 __m128d fix, fiy, fiz;
781 __m128d rbai, rbaj, f_gb, f_gb_ai;
782 __m128d xmm1, xmm2, xmm3;
784 const __m128d two = _mm_set1_pd(2.0);
790 /* Loop to get the proper form for the Born radius term, sse style */
794 if (gb_algorithm == egbSTILL)
796 for (i = n0; i < n1; i++)
799 rb[i] = (2 * rbi * rbi * dvda[i])/ONE_4PI_EPS0;
802 else if (gb_algorithm == egbHCT)
804 for (i = n0; i < n1; i++)
807 rb[i] = rbi * rbi * dvda[i];
810 else if (gb_algorithm == egbOBC)
812 for (i = n0; i < n1; i++)
815 rb[i] = rbi * rbi * born->drobc[i] * dvda[i];
819 jz = _mm_setzero_pd();
823 for (i = 0; i < nl->nri; i++)
827 is3 = 3*nl->shift[i];
829 shY = shiftvec[is3+1];
830 shZ = shiftvec[is3+2];
832 nj1 = nl->jindex[i+1];
834 ix = _mm_set1_pd(shX+x[ii3+0]);
835 iy = _mm_set1_pd(shY+x[ii3+1]);
836 iz = _mm_set1_pd(shZ+x[ii3+2]);
838 rbai = _mm_load1_pd(rb+ii);
839 fix = _mm_setzero_pd();
840 fiy = _mm_setzero_pd();
841 fiz = _mm_setzero_pd();
844 for (k = nj0; k < nj1-1; k += 2)
852 GMX_MM_LOAD_1RVEC_2POINTERS_PD(x+j3A, x+j3B, jx, jy, jz);
854 dx = _mm_sub_pd(ix, jx);
855 dy = _mm_sub_pd(iy, jy);
856 dz = _mm_sub_pd(iz, jz);
858 GMX_MM_LOAD_2VALUES_PD(rb+jnrA, rb+jnrB, rbaj);
860 /* load chain rule terms for j1-4 */
861 f_gb = _mm_load_pd(dadx);
863 f_gb_ai = _mm_load_pd(dadx);
866 /* calculate scalar force */
867 f_gb = _mm_mul_pd(f_gb, rbai);
868 f_gb_ai = _mm_mul_pd(f_gb_ai, rbaj);
869 f_gb = _mm_add_pd(f_gb, f_gb_ai);
871 tx = _mm_mul_pd(f_gb, dx);
872 ty = _mm_mul_pd(f_gb, dy);
873 tz = _mm_mul_pd(f_gb, dz);
875 fix = _mm_add_pd(fix, tx);
876 fiy = _mm_add_pd(fiy, ty);
877 fiz = _mm_add_pd(fiz, tz);
879 GMX_MM_DECREMENT_1RVEC_2POINTERS_PD(f+j3A, f+j3B, tx, ty, tz);
882 /*deal with odd elements */
888 GMX_MM_LOAD_1RVEC_1POINTER_PD(x+j3A, jx, jy, jz);
890 dx = _mm_sub_sd(ix, jx);
891 dy = _mm_sub_sd(iy, jy);
892 dz = _mm_sub_sd(iz, jz);
894 GMX_MM_LOAD_1VALUE_PD(rb+jnrA, rbaj);
896 /* load chain rule terms */
897 f_gb = _mm_load_pd(dadx);
899 f_gb_ai = _mm_load_pd(dadx);
902 /* calculate scalar force */
903 f_gb = _mm_mul_sd(f_gb, rbai);
904 f_gb_ai = _mm_mul_sd(f_gb_ai, rbaj);
905 f_gb = _mm_add_sd(f_gb, f_gb_ai);
907 tx = _mm_mul_sd(f_gb, dx);
908 ty = _mm_mul_sd(f_gb, dy);
909 tz = _mm_mul_sd(f_gb, dz);
911 fix = _mm_add_sd(fix, tx);
912 fiy = _mm_add_sd(fiy, ty);
913 fiz = _mm_add_sd(fiz, tz);
915 GMX_MM_DECREMENT_1RVEC_1POINTER_PD(f+j3A, tx, ty, tz);
918 /* fix/fiy/fiz now contain four partial force terms, that all should be
919 * added to the i particle forces and shift forces.
921 gmx_mm_update_iforce_1atom_pd(&fix, &fiy, &fiz, f+ii3, fshift+is3);
928 /* keep compiler happy */
929 int genborn_sse2_dummy;
931 #endif /* SSE2 intrinsics available */