?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Updated intel syntax x86-64 asm files to also support MS win64 call convention (ifdef...
[gromacs/rigid-bodies.git] / src / gmxlib / nonbonded / nonbonded.c
blob08f88f62807264e1973c3c8eacb215f26129cc85
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * GROningen Mixture of Alchemy and Childrens' Stories
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #ifdef GMX_THREAD_SHM_FDECOMP
41 #include <thread_mpi.h>
42 #endif
43
44
45 #include <stdio.h>
46 #include <stdlib.h>
47 #include "typedefs.h"
48 #include "txtdump.h"
49 #include "smalloc.h"
50 #include "ns.h"
51 #include "vec.h"
52 #include "maths.h"
53 #include "macros.h"
54 #include "force.h"
55 #include "names.h"
56 #include "main.h"
57 #include "xvgr.h"
58 #include "gmx_fatal.h"
59 #include "physics.h"
60 #include "force.h"
61 #include "bondf.h"
62 #include "nrnb.h"
63 #include "smalloc.h"
64 #include "nonbonded.h"
65
66 #include "nb_kernel_c/nb_kernel_c.h"
67 #include "nb_free_energy.h"
68 #include "nb_generic.h"
69 #include "nb_generic_cg.h"
70
71
72 /* 1,4 interactions uses kernel 330 directly */
73 #include "nb_kernel_c/nb_kernel330.h"
74
75 #ifdef GMX_PPC_ALTIVEC
76 #include "nb_kernel_ppc_altivec/nb_kernel_ppc_altivec.h"
77 #endif
78
79 #if defined(GMX_IA32_SSE)
80 #include "nb_kernel_ia32_sse/nb_kernel_ia32_sse.h"
81 #endif
82
83 #if defined(GMX_IA32_SSE2)
84 #include "nb_kernel_ia32_sse2/nb_kernel_ia32_sse2.h"
85 #endif
86
87 #if defined(GMX_X86_64_SSE)
88 #include "nb_kernel_x86_64_sse/nb_kernel_x86_64_sse.h"
89 #endif
90
91 #if defined(GMX_X86_64_SSE2)
92 #include "nb_kernel_x86_64_sse2/nb_kernel_x86_64_sse2.h"
93 #endif
94
95 #if defined(GMX_SSE2)
96 # ifdef GMX_DOUBLE
97 # include "nb_kernel_sse2_double/nb_kernel_sse2_double.h"
98 # else
99 # include "nb_kernel_sse2_single/nb_kernel_sse2_single.h"
100 # endif
101 #endif
102
103 #if defined(GMX_FORTRAN)
104 # ifdef GMX_DOUBLE
105 # include "nb_kernel_f77_double/nb_kernel_f77_double.h"
106 # else
107 # include "nb_kernel_f77_single/nb_kernel_f77_single.h"
108 # endif
109 #endif
110
111 #if (defined GMX_IA64_ASM && defined GMX_DOUBLE)
112 #include "nb_kernel_ia64_double/nb_kernel_ia64_double.h"
113 #endif
114
115 #if (defined GMX_IA64_ASM && !defined GMX_DOUBLE)
116 #include "nb_kernel_ia64_single/nb_kernel_ia64_single.h"
117 #endif
118
119 #ifdef GMX_POWER6
120 #include "nb_kernel_power6/nb_kernel_power6.h"
121 #endif
122
123 #ifdef GMX_BLUEGENE
124 #include "nb_kernel_bluegene/nb_kernel_bluegene.h"
125 #endif
126
127
128
129 enum { TABLE_NONE, TABLE_COMBINED, TABLE_COUL, TABLE_VDW, TABLE_NR };
130
131
132 /* Table version for each kernel.
133 */
134 static const int
135 nb_kernel_table[eNR_NBKERNEL_NR] =
136 {
137 TABLE_NONE, /* kernel010 */
138 TABLE_NONE, /* kernel020 */
139 TABLE_VDW, /* kernel030 */
140 TABLE_NONE, /* kernel100 */
141 TABLE_NONE, /* kernel101 */
142 TABLE_NONE, /* kernel102 */
143 TABLE_NONE, /* kernel103 */
144 TABLE_NONE, /* kernel104 */
145 TABLE_NONE, /* kernel110 */
146 TABLE_NONE, /* kernel111 */
147 TABLE_NONE, /* kernel112 */
148 TABLE_NONE, /* kernel113 */
149 TABLE_NONE, /* kernel114 */
150 TABLE_NONE, /* kernel120 */
151 TABLE_NONE, /* kernel121 */
152 TABLE_NONE, /* kernel122 */
153 TABLE_NONE, /* kernel123 */
154 TABLE_NONE, /* kernel124 */
155 TABLE_VDW, /* kernel130 */
156 TABLE_VDW, /* kernel131 */
157 TABLE_VDW, /* kernel132 */
158 TABLE_VDW, /* kernel133 */
159 TABLE_VDW, /* kernel134 */
160 TABLE_NONE, /* kernel200 */
161 TABLE_NONE, /* kernel201 */
162 TABLE_NONE, /* kernel202 */
163 TABLE_NONE, /* kernel203 */
164 TABLE_NONE, /* kernel204 */
165 TABLE_NONE, /* kernel210 */
166 TABLE_NONE, /* kernel211 */
167 TABLE_NONE, /* kernel212 */
168 TABLE_NONE, /* kernel213 */
169 TABLE_NONE, /* kernel214 */
170 TABLE_NONE, /* kernel220 */
171 TABLE_NONE, /* kernel221 */
172 TABLE_NONE, /* kernel222 */
173 TABLE_NONE, /* kernel223 */
174 TABLE_NONE, /* kernel224 */
175 TABLE_VDW, /* kernel230 */
176 TABLE_VDW, /* kernel231 */
177 TABLE_VDW, /* kernel232 */
178 TABLE_VDW, /* kernel233 */
179 TABLE_VDW, /* kernel234 */
180 TABLE_COUL, /* kernel300 */
181 TABLE_COUL, /* kernel301 */
182 TABLE_COUL, /* kernel302 */
183 TABLE_COUL, /* kernel303 */
184 TABLE_COUL, /* kernel304 */
185 TABLE_COUL, /* kernel310 */
186 TABLE_COUL, /* kernel311 */
187 TABLE_COUL, /* kernel312 */
188 TABLE_COUL, /* kernel313 */
189 TABLE_COUL, /* kernel314 */
190 TABLE_COUL, /* kernel320 */
191 TABLE_COUL, /* kernel321 */
192 TABLE_COUL, /* kernel322 */
193 TABLE_COUL, /* kernel323 */
194 TABLE_COUL, /* kernel324 */
195 TABLE_COMBINED, /* kernel330 */
196 TABLE_COMBINED, /* kernel331 */
197 TABLE_COMBINED, /* kernel332 */
198 TABLE_COMBINED, /* kernel333 */
199 TABLE_COMBINED, /* kernel334 */
200 TABLE_NONE, /* kernel400 */
201 TABLE_NONE, /* kernel410 */
202 TABLE_VDW /* kernel430 */
203 };
204
205
206
207 static nb_kernel_t **
208 nb_kernel_list = NULL;
209
210
211 void
212 gmx_setup_kernels(FILE *fplog,bool bGenericKernelOnly)
213 {
214 int i;
215
216 snew(nb_kernel_list,eNR_NBKERNEL_NR);
217
218 /* Note that later calls overwrite earlier, so the preferred (fastest)
219 * version should be at the end. For instance, we call SSE after 3DNow.
220 */
221
222 for(i=0; i<eNR_NBKERNEL_NR; i++)
223 {
224 nb_kernel_list[i] = NULL;
225 }
226
227 if (bGenericKernelOnly)
228 {
229 return;
230 }
231
232 if(fplog)
233 {
234 fprintf(fplog,"Configuring nonbonded kernels...\n");
235 }
236
237 nb_kernel_setup(fplog,nb_kernel_list);
238
239 if(getenv("GMX_NOOPTIMIZEDKERNELS") != NULL)
240 {
241 return;
242 }
243
244 /* Setup kernels. The last called setup routine will overwrite earlier assignments,
245 * so we should e.g. test SSE3 support _after_ SSE2 support,
246 * and call e.g. Fortran setup before SSE.
247 */
248
249 #if defined(GMX_FORTRAN) && defined(GMX_DOUBLE)
250 nb_kernel_setup_f77_double(fplog,nb_kernel_list);
251 #endif
252
253 #if defined(GMX_FORTRAN) && !defined(GMX_DOUBLE)
254 nb_kernel_setup_f77_single(fplog,nb_kernel_list);
255 #endif
256
257 #ifdef GMX_BLUEGENE
258 nb_kernel_setup_bluegene(fplog,nb_kernel_list);
259 #endif
260
261 #ifdef GMX_POWER6
262 nb_kernel_setup_power6(fplog,nb_kernel_list);
263 #endif
264
265 #ifdef GMX_PPC_ALTIVEC
266 nb_kernel_setup_ppc_altivec(fplog,nb_kernel_list);
267 #endif
268
269 #if defined(GMX_IA32_SSE)
270 nb_kernel_setup_ia32_sse(fplog,nb_kernel_list);
271 #endif
272
273 #if defined(GMX_IA32_SSE2)
274 nb_kernel_setup_ia32_sse2(fplog,nb_kernel_list);
275 #endif
276
277 #if defined(GMX_X86_64_SSE)
278 nb_kernel_setup_x86_64_sse(fplog,nb_kernel_list);
279 #endif
280
281 #if defined(GMX_X86_64_SSE2)
282 nb_kernel_setup_x86_64_sse2(fplog,nb_kernel_list);
283 #endif
284
285 #if (defined GMX_IA64_ASM && defined GMX_DOUBLE)
286 nb_kernel_setup_ia64_double(fplog,nb_kernel_list);
287 #endif
288
289 #if (defined GMX_IA64_ASM && !defined GMX_DOUBLE)
290 nb_kernel_setup_ia64_single(fplog,nb_kernel_list);
291 #endif
292
293 if(fplog)
294 {
295 fprintf(fplog,"\n\n");
296 }
297 }
298
299
300 void do_nonbonded(t_commrec *cr,t_forcerec *fr,
301 rvec x[],rvec f[],t_mdatoms *mdatoms,t_blocka *excl,
302 real egnb[],real egcoul[],real egpol[],rvec box_size,
303 t_nrnb *nrnb,real lambda,real *dvdlambda,
304 int nls,int eNL,int flags)
305 {
306 bool bLR,bDoForces,bForeignLambda;
307 t_nblist * nlist;
308 real * fshift;
309 int n,n0,n1,i,i0,i1,nrnb_ind,sz;
310 t_nblists *nblists;
311 bool bWater;
312 nb_kernel_t * kernelptr;
313 FILE * fp;
314 int fac=0;
315 int nthreads = 1;
316 int tabletype;
317 int outeriter,inneriter;
318 real * tabledata = NULL;
319 gmx_gbdata_t gbdata;
320
321 bLR = (flags & GMX_DONB_LR);
322 bDoForces = (flags & GMX_DONB_FORCES);
323 bForeignLambda = (flags & GMX_DONB_FOREIGNLAMBDA);
324
325 gbdata.gb_epsilon_solvent = fr->gb_epsilon_solvent;
326 gbdata.epsilon_r = fr->epsilon_r;
327 gbdata.gpol = egpol;
328
329 if(fr->bAllvsAll)
330 {
331 if(fr->bGB)
332 {
333 #if (defined GMX_SSE2 || defined GMX_X86_64_SSE || defined GMX_X86_64_SSE2 || defined GMX_IA32_SSE || defined GMX_IA32_SSE2)
334 # ifdef GMX_DOUBLE
335 if(fr->UseOptimizedKernels)
336 {
337 nb_kernel_allvsallgb_sse2_double(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
338 &outeriter,&inneriter,&fr->AllvsAll_work);
339 }
340 else
341 {
342 nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
343 &outeriter,&inneriter,&fr->AllvsAll_work);
344 }
345 # else /* not double */
346 if(fr->UseOptimizedKernels)
347 {
348 nb_kernel_allvsallgb_sse2_single(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
349 &outeriter,&inneriter,&fr->AllvsAll_work);
350 }
351 else
352 {
353 nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
354 &outeriter,&inneriter,&fr->AllvsAll_work);
355 }
356 # endif /* double/single alt. */
357 #else /* no SSE support compiled in */
358 nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
359 &outeriter,&inneriter,&fr->AllvsAll_work);
360 #endif
361 inc_nrnb(nrnb,eNR_NBKERNEL_ALLVSALLGB,inneriter);
362 }
363 else
364 {
365 #if (defined GMX_SSE2 || defined GMX_X86_64_SSE || defined GMX_X86_64_SSE2 || defined GMX_IA32_SSE || defined GMX_IA32_SSE2)
366 # ifdef GMX_DOUBLE
367 if(fr->UseOptimizedKernels)
368 {
369 nb_kernel_allvsall_sse2_double(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
370 &outeriter,&inneriter,&fr->AllvsAll_work);
371 }
372 else
373 {
374 nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
375 &outeriter,&inneriter,&fr->AllvsAll_work);
376 }
377
378 # else /* not double */
379 if(fr->UseOptimizedKernels)
380 {
381 nb_kernel_allvsall_sse2_single(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
382 &outeriter,&inneriter,&fr->AllvsAll_work);
383 }
384 else
385 {
386 nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
387 &outeriter,&inneriter,&fr->AllvsAll_work);
388 }
389
390 # endif /* double/single check */
391 #else /* No SSE2 support compiled in */
392 nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
393 &outeriter,&inneriter,&fr->AllvsAll_work);
394 #endif
395
396 inc_nrnb(nrnb,eNR_NBKERNEL_ALLVSALL,inneriter);
397 }
398 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,outeriter);
399 return;
400 }
401
402 if (eNL >= 0)
403 {
404 i0 = eNL;
405 i1 = i0+1;
406 }
407 else
408 {
409 i0 = 0;
410 i1 = eNL_NR;
411 }
412
413 if (nls >= 0)
414 {
415 n0 = nls;
416 n1 = nls+1;
417 }
418 else
419 {
420 n0 = 0;
421 n1 = fr->nnblists;
422 }
423
424 if(nb_kernel_list == NULL)
425 {
426 gmx_fatal(FARGS,"gmx_setup_kernels has not been called");
427 }
428
429 fshift = fr->fshift[0];
430
431 for(n=n0; (n<n1); n++)
432 {
433 nblists = &fr->nblists[n];
434 for(i=i0; (i<i1); i++)
435 {
436 outeriter = inneriter = 0;
437
438 if (bLR)
439 {
440 nlist = &(nblists->nlist_lr[i]);
441 }
442 else
443 {
444 nlist = &(nblists->nlist_sr[i]);
445 }
446
447 if (nlist->nri > 0)
448 {
449 nrnb_ind = nlist->il_code;
450
451 if(nrnb_ind==eNR_NBKERNEL_FREE_ENERGY)
452 {
453 /* generic free energy, use combined table */
454 tabledata = nblists->tab.tab;
455 }
456 else
457 {
458 if (bForeignLambda)
459 {
460 /* We don't need the non-perturbed interactions */
461 continue;
462 }
463
464 tabletype = nb_kernel_table[nrnb_ind];
465
466 /* normal kernels, not free energy */
467 if (!bDoForces)
468 {
469 nrnb_ind += eNR_NBKERNEL_NR/2;
470 }
471
472 if(tabletype == TABLE_COMBINED)
473 {
474 tabledata = nblists->tab.tab;
475 }
476 else if(tabletype == TABLE_COUL)
477 {
478 tabledata = nblists->coultab;
479 }
480 else if(tabletype == TABLE_VDW)
481 {
482 tabledata = nblists->vdwtab;
483 }
484 else
485 {
486 tabledata = NULL;
487 }
488 }
489
490 nlist->count = 0;
491
492 if(nlist->free_energy)
493 {
494 if(nlist->ivdw==2)
495 {
496 gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
497 }
498
499 gmx_nb_free_energy_kernel(nlist->icoul,
500 nlist->ivdw,
501 nlist->nri,
502 nlist->iinr,
503 nlist->jindex,
504 nlist->jjnr,
505 nlist->shift,
506 fr->shift_vec[0],
507 fshift,
508 nlist->gid,
509 x[0],
510 f[0],
511 mdatoms->chargeA,
512 mdatoms->chargeB,
513 fr->epsfac,
514 fr->k_rf,
515 fr->c_rf,
516 fr->ewaldcoeff,
517 egcoul,
518 mdatoms->typeA,
519 mdatoms->typeB,
520 fr->ntype,
521 fr->nbfp,
522 egnb,
523 nblists->tab.scale,
524 tabledata,
525 lambda,
526 dvdlambda,
527 fr->sc_alpha,
528 fr->sc_power,
529 fr->sc_sigma6_def,
530 fr->sc_sigma6_min,
531 bDoForces,
532 &outeriter,
533 &inneriter);
534 }
535 else if (nlist->enlist == enlistCG_CG)
536 {
537 /* Call the charge group based inner loop */
538 gmx_nb_generic_cg_kernel(nlist,
539 fr,
540 mdatoms,
541 x[0],
542 f[0],
543 fshift,
544 egcoul,
545 egnb,
546 nblists->tab.scale,
547 tabledata,
548 &outeriter,
549 &inneriter);
550 }
551 else
552 {
553 /* Not free energy */
554
555 kernelptr = nb_kernel_list[nrnb_ind];
556
557 if (kernelptr == NULL)
558 {
559 /* Call a generic nonbonded kernel */
560
561 /* If you want to hack/test your own interactions,
562 * do it in this routine and make sure it is called
563 * by setting the environment variable GMX_NB_GENERIC.
564 */
565 gmx_nb_generic_kernel(nlist,
566 fr,
567 mdatoms,
568 x[0],
569 f[0],
570 fshift,
571 egcoul,
572 egnb,
573 nblists->tab.scale,
574 tabledata,
575 &outeriter,
576 &inneriter);
577 }
578 else
579 {
580 /* Call nonbonded kernel from function pointer */
581
582 (*kernelptr)( &(nlist->nri),
583 nlist->iinr,
584 nlist->jindex,
585 nlist->jjnr,
586 nlist->shift,
587 fr->shift_vec[0],
588 fshift,
589 nlist->gid,
590 x[0],
591 f[0],
592 mdatoms->chargeA,
593 &(fr->epsfac),
594 &(fr->k_rf),
595 &(fr->c_rf),
596 egcoul,
597 mdatoms->typeA,
598 &(fr->ntype),
599 fr->nbfp,
600 egnb,
601 &(nblists->tab.scale),
602 tabledata,
603 fr->invsqrta,
604 fr->dvda,
605 &(fr->gbtabscale),
606 fr->gbtab.tab,
607 &nthreads,
608 &(nlist->count),
609 nlist->mtx,
610 &outeriter,
611 &inneriter,
612 (real *)&gbdata);
613 }
614 }
615
616 /* Update flop accounting */
617
618 /* Outer loop in kernel */
619 switch (nlist->enlist) {
620 case enlistATOM_ATOM: fac = 1; break;
621 case enlistSPC_ATOM: fac = 3; break;
622 case enlistSPC_SPC: fac = 9; break;
623 case enlistTIP4P_ATOM: fac = 4; break;
624 case enlistTIP4P_TIP4P: fac = 16; break;
625 case enlistCG_CG: fac = 1; break;
626 }
627 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,fac*outeriter);
628
629 /* inner loop in kernel */
630 inc_nrnb(nrnb,nrnb_ind,inneriter);
631 }
632 }
633 }
634 }
635
636
637 real
638 do_listed_vdw_q(int ftype,int nbonds,
639 const t_iatom iatoms[],const t_iparams iparams[],
640 const rvec x[],rvec f[],rvec fshift[],
641 const t_pbc *pbc,const t_graph *g,
642 real lambda,real *dvdlambda,
643 const t_mdatoms *md,
644 const t_forcerec *fr,gmx_grppairener_t *grppener,
645 int *global_atom_index)
646 {
647 static bool bWarn=FALSE;
648 real eps,r2,*tab,rtab2=0;
649 rvec dx,x14[2],f14[2];
650 int i,ai,aj,itype;
651 int typeA[2]={0,0},typeB[2]={0,1};
652 real chargeA[2]={0,0},chargeB[2];
653 int gid,shift_vir,shift_f;
654 int j_index[] = { 0, 1 };
655 int i0=0,i1=1,i2=2;
656 ivec dt;
657 int outeriter,inneriter;
658 int nthreads = 1;
659 int count;
660 real krf,crf,tabscale;
661 int ntype=0;
662 real *nbfp=NULL;
663 real *egnb=NULL,*egcoul=NULL;
664 t_nblist tmplist;
665 int icoul,ivdw;
666 bool bMolPBC,bFreeEnergy;
667
668 #if GMX_THREAD_SHM_FDECOMP
669 pthread_mutex_t mtx;
670 #else
671 void * mtx = NULL;
672 #endif
673
674
675 #if GMX_THREAD_SHM_FDECOMP
676 pthread_mutex_initialize(&mtx);
677 #endif
678
679 bMolPBC = fr->bMolPBC;
680
681 switch (ftype) {
682 case F_LJ14:
683 case F_LJC14_Q:
684 eps = fr->epsfac*fr->fudgeQQ;
685 ntype = 1;
686 egnb = grppener->ener[egLJ14];
687 egcoul = grppener->ener[egCOUL14];
688 break;
689 case F_LJC_PAIRS_NB:
690 eps = fr->epsfac;
691 ntype = 1;
692 egnb = grppener->ener[egLJSR];
693 egcoul = grppener->ener[egCOULSR];
694 break;
695 default:
696 gmx_fatal(FARGS,"Unknown function type %d in do_nonbonded14",
697 ftype);
698 }
699 tab = fr->tab14.tab;
700 rtab2 = sqr(fr->tab14.r);
701 tabscale = fr->tab14.scale;
702
703 krf = fr->k_rf;
704 crf = fr->c_rf;
705
706 /* Determine the values for icoul/ivdw. */
707 if (fr->bEwald) {
708 icoul = 1;
709 }
710 else if(fr->bcoultab)
711 {
712 icoul = 3;
713 }
714 else if(fr->eeltype == eelRF_NEC)
715 {
716 icoul = 2;
717 }
718 else
719 {
720 icoul = 1;
721 }
722
723 if(fr->bvdwtab)
724 {
725 ivdw = 3;
726 }
727 else if(fr->bBHAM)
728 {
729 ivdw = 2;
730 }
731 else
732 {
733 ivdw = 1;
734 }
735
736
737 /* We don't do SSE or altivec here, due to large overhead for 4-fold
738 * unrolling on short lists
739 */
740
741 bFreeEnergy = FALSE;
742 for(i=0; (i<nbonds); )
743 {
744 itype = iatoms[i++];
745 ai = iatoms[i++];
746 aj = iatoms[i++];
747 gid = GID(md->cENER[ai],md->cENER[aj],md->nenergrp);
748
749 switch (ftype) {
750 case F_LJ14:
751 bFreeEnergy =
752 (fr->efep != efepNO &&
753 ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
754 iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
755 iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
756 chargeA[0] = md->chargeA[ai];
757 chargeA[1] = md->chargeA[aj];
758 nbfp = (real *)&(iparams[itype].lj14.c6A);
759 break;
760 case F_LJC14_Q:
761 eps = fr->epsfac*iparams[itype].ljc14.fqq;
762 chargeA[0] = iparams[itype].ljc14.qi;
763 chargeA[1] = iparams[itype].ljc14.qj;
764 nbfp = (real *)&(iparams[itype].ljc14.c6);
765 break;
766 case F_LJC_PAIRS_NB:
767 chargeA[0] = iparams[itype].ljcnb.qi;
768 chargeA[1] = iparams[itype].ljcnb.qj;
769 nbfp = (real *)&(iparams[itype].ljcnb.c6);
770 break;
771 }
772
773 if (!bMolPBC)
774 {
775 /* This is a bonded interaction, atoms are in the same box */
776 shift_f = CENTRAL;
777 r2 = distance2(x[ai],x[aj]);
778 }
779 else
780 {
781 /* Apply full periodic boundary conditions */
782 shift_f = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
783 r2 = norm2(dx);
784 }
785
786 if (r2 >= rtab2)
787 {
788 if (!bWarn)
789 {
790 fprintf(stderr,"Warning: 1-4 interaction between %d and %d "
791 "at distance %.3f which is larger than the 1-4 table size %.3f nm\n",
792 glatnr(global_atom_index,ai),
793 glatnr(global_atom_index,aj),
794 sqrt(r2), sqrt(rtab2));
795 fprintf(stderr,"These are ignored for the rest of the simulation\n");
796 fprintf(stderr,"This usually means your system is exploding,\n"
797 "if not, you should increase table-extension in your mdp file\n"
798 "or with user tables increase the table size\n");
799 bWarn = TRUE;
800 }
801 if (debug)
802 fprintf(debug,"%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
803 x[ai][XX],x[ai][YY],x[ai][ZZ],
804 x[aj][XX],x[aj][YY],x[aj][ZZ],
805 glatnr(global_atom_index,ai),
806 glatnr(global_atom_index,aj),
807 sqrt(r2));
808 }
809 else
810 {
811 copy_rvec(x[ai],x14[0]);
812 copy_rvec(x[aj],x14[1]);
813 clear_rvec(f14[0]);
814 clear_rvec(f14[1]);
815 #ifdef DEBUG
816 fprintf(debug,"LJ14: grp-i=%2d, grp-j=%2d, ngrp=%2d, GID=%d\n",
817 md->cENER[ai],md->cENER[aj],md->nenergrp,gid);
818 #endif
819
820 outeriter = inneriter = count = 0;
821 if (bFreeEnergy)
822 {
823 chargeB[0] = md->chargeB[ai];
824 chargeB[1] = md->chargeB[aj];
825 /* We pass &(iparams[itype].lj14.c6A) as LJ parameter matrix
826 * to the innerloops.
827 * Here we use that the LJ-14 parameters are stored in iparams
828 * as c6A,c12A,c6B,c12B, which are referenced correctly
829 * in the innerloops if we assign type combinations 0-0 and 0-1
830 * to atom pair ai-aj in topologies A and B respectively.
831 */
832 if(ivdw==2)
833 {
834 gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
835 }
836 count = 0;
837 gmx_nb_free_energy_kernel(icoul,
838 ivdw,
839 i1,
840 &i0,
841 j_index,
842 &i1,
843 &shift_f,
844 fr->shift_vec[0],
845 fshift[0],
846 &gid,
847 x14[0],
848 f14[0],
849 chargeA,
850 chargeB,
851 eps,
852 krf,
853 crf,
854 fr->ewaldcoeff,
855 egcoul,
856 typeA,
857 typeB,
858 ntype,
859 nbfp,
860 egnb,
861 tabscale,
862 tab,
863 lambda,
864 dvdlambda,
865 fr->sc_alpha,
866 fr->sc_power,
867 fr->sc_sigma6_def,
868 fr->sc_sigma6_min,
869 TRUE,
870 &outeriter,
871 &inneriter);
872 }
873 else
874 {
875 /* Not perturbed - call kernel 330 */
876 nb_kernel330
877 ( &i1,
878 &i0,
879 j_index,
880 &i1,
881 &shift_f,
882 fr->shift_vec[0],
883 fshift[0],
884 &gid,
885 x14[0],
886 f14[0],
887 chargeA,
888 &eps,
889 &krf,
890 &crf,
891 egcoul,
892 typeA,
893 &ntype,
894 nbfp,
895 egnb,
896 &tabscale,
897 tab,
898 NULL,
899 NULL,
900 NULL,
901 NULL,
902 &nthreads,
903 &count,
904 (void *)&mtx,
905 &outeriter,
906 &inneriter,
907 NULL);
908 }
909
910 /* Add the forces */
911 rvec_inc(f[ai],f14[0]);
912 rvec_dec(f[aj],f14[0]);
913
914 if (g)
915 {
916 /* Correct the shift forces using the graph */
917 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
918 shift_vir = IVEC2IS(dt);
919 rvec_inc(fshift[shift_vir],f14[0]);
920 rvec_dec(fshift[CENTRAL],f14[0]);
921 }
922
923 /* flops: eNR_KERNEL_OUTER + eNR_KERNEL330 + 12 */
924 }
925 }
926 return 0.0;
927 }
928
929