?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
added Verlet scheme and NxN non-bonded functionality
[gromacs.git] / src / gmxlib / nonbonded / nonbonded.c
blobd0120c03fa88385b586da6327ca2be438573df8e
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_kernel_adress_c/nb_kernel_c_adress.h"
68 #include "nb_free_energy.h"
69 #include "nb_generic.h"
70 #include "nb_generic_cg.h"
71 #include "nb_generic_adress.h"
72
73
74 /* 1,4 interactions uses kernel 330 directly */
75 #include "nb_kernel_c/nb_kernel330.h"
76 #include "nb_kernel_adress_c/nb_kernel330_adress.h"
77
78
79
80 #if 0 && defined (GMX_X86_SSE2)
81 # ifdef GMX_DOUBLE
82 # include "nb_kernel_sse2_double/nb_kernel_sse2_double.h"
83 # else
84 # include "nb_kernel_sse2_single/nb_kernel_sse2_single.h"
85 # endif
86 #endif
87
88 #if defined(GMX_FORTRAN)
89 # ifdef GMX_DOUBLE
90 # include "nb_kernel_f77_double/nb_kernel_f77_double.h"
91 # else
92 # include "nb_kernel_f77_single/nb_kernel_f77_single.h"
93 # endif
94 #endif
95
96
97 #ifdef GMX_POWER6
98 #include "nb_kernel_power6/nb_kernel_power6.h"
99 #endif
100
101 #ifdef GMX_BLUEGENE
102 #include "nb_kernel_bluegene/nb_kernel_bluegene.h"
103 #endif
104
105
106
107 enum { TABLE_NONE, TABLE_COMBINED, TABLE_COUL, TABLE_VDW, TABLE_NR };
108
109
110 /* Table version for each kernel.
111 */
112 static const int
113 nb_kernel_table[eNR_NBKERNEL_NR] =
114 {
115 TABLE_NONE, /* kernel010 */
116 TABLE_NONE, /* kernel020 */
117 TABLE_VDW, /* kernel030 */
118 TABLE_NONE, /* kernel100 */
119 TABLE_NONE, /* kernel101 */
120 TABLE_NONE, /* kernel102 */
121 TABLE_NONE, /* kernel103 */
122 TABLE_NONE, /* kernel104 */
123 TABLE_NONE, /* kernel110 */
124 TABLE_NONE, /* kernel111 */
125 TABLE_NONE, /* kernel112 */
126 TABLE_NONE, /* kernel113 */
127 TABLE_NONE, /* kernel114 */
128 TABLE_NONE, /* kernel120 */
129 TABLE_NONE, /* kernel121 */
130 TABLE_NONE, /* kernel122 */
131 TABLE_NONE, /* kernel123 */
132 TABLE_NONE, /* kernel124 */
133 TABLE_VDW, /* kernel130 */
134 TABLE_VDW, /* kernel131 */
135 TABLE_VDW, /* kernel132 */
136 TABLE_VDW, /* kernel133 */
137 TABLE_VDW, /* kernel134 */
138 TABLE_NONE, /* kernel200 */
139 TABLE_NONE, /* kernel201 */
140 TABLE_NONE, /* kernel202 */
141 TABLE_NONE, /* kernel203 */
142 TABLE_NONE, /* kernel204 */
143 TABLE_NONE, /* kernel210 */
144 TABLE_NONE, /* kernel211 */
145 TABLE_NONE, /* kernel212 */
146 TABLE_NONE, /* kernel213 */
147 TABLE_NONE, /* kernel214 */
148 TABLE_NONE, /* kernel220 */
149 TABLE_NONE, /* kernel221 */
150 TABLE_NONE, /* kernel222 */
151 TABLE_NONE, /* kernel223 */
152 TABLE_NONE, /* kernel224 */
153 TABLE_VDW, /* kernel230 */
154 TABLE_VDW, /* kernel231 */
155 TABLE_VDW, /* kernel232 */
156 TABLE_VDW, /* kernel233 */
157 TABLE_VDW, /* kernel234 */
158 TABLE_COUL, /* kernel300 */
159 TABLE_COUL, /* kernel301 */
160 TABLE_COUL, /* kernel302 */
161 TABLE_COUL, /* kernel303 */
162 TABLE_COUL, /* kernel304 */
163 TABLE_COUL, /* kernel310 */
164 TABLE_COUL, /* kernel311 */
165 TABLE_COUL, /* kernel312 */
166 TABLE_COUL, /* kernel313 */
167 TABLE_COUL, /* kernel314 */
168 TABLE_COUL, /* kernel320 */
169 TABLE_COUL, /* kernel321 */
170 TABLE_COUL, /* kernel322 */
171 TABLE_COUL, /* kernel323 */
172 TABLE_COUL, /* kernel324 */
173 TABLE_COMBINED, /* kernel330 */
174 TABLE_COMBINED, /* kernel331 */
175 TABLE_COMBINED, /* kernel332 */
176 TABLE_COMBINED, /* kernel333 */
177 TABLE_COMBINED, /* kernel334 */
178 TABLE_NONE, /* kernel400 */
179 TABLE_NONE, /* kernel410 */
180 TABLE_VDW /* kernel430 */
181 };
182
183
184
185 static nb_kernel_t **
186 nb_kernel_list = NULL;
187
188 static nb_adress_kernel_t **
189 nb_kernel_list_adress = NULL;
190
191 void
192 gmx_setup_kernels(FILE *fplog,t_forcerec *fr,gmx_bool bGenericKernelOnly)
193 {
194 int i;
195
196 snew(nb_kernel_list,eNR_NBKERNEL_NR);
197
198 /* Note that later calls overwrite earlier, so the preferred (fastest)
199 * version should be at the end. For instance, we call SSE after 3DNow.
200 */
201
202 for(i=0; i<eNR_NBKERNEL_NR; i++)
203 {
204 nb_kernel_list[i] = NULL;
205 }
206
207 if (bGenericKernelOnly)
208 {
209 return;
210 }
211
212 if(fplog)
213 {
214 fprintf(fplog,"Configuring nonbonded kernels...\n");
215 }
216
217 nb_kernel_setup(fplog,nb_kernel_list);
218
219 if(fr->use_cpu_acceleration==FALSE)
220 {
221 return;
222 }
223
224 /* Setup kernels. The last called setup routine will overwrite earlier assignments,
225 * so we should e.g. test SSE3 support _after_ SSE2 support,
226 * and call e.g. Fortran setup before SSE.
227 */
228
229 #if defined(GMX_FORTRAN) && defined(GMX_DOUBLE)
230 nb_kernel_setup_f77_double(fplog,nb_kernel_list);
231 #endif
232
233 #if defined(GMX_FORTRAN) && !defined(GMX_DOUBLE)
234 nb_kernel_setup_f77_single(fplog,nb_kernel_list);
235 #endif
236
237 #ifdef GMX_BLUEGENE
238 nb_kernel_setup_bluegene(fplog,nb_kernel_list);
239 #endif
240
241 #ifdef GMX_POWER6
242 nb_kernel_setup_power6(fplog,nb_kernel_list);
243 #endif
244
245 if(fplog)
246 {
247 fprintf(fplog,"\n\n");
248 }
249 }
250
251 void
252 gmx_setup_adress_kernels(FILE *fplog,gmx_bool bGenericKernelOnly) {
253 int i;
254
255 snew(nb_kernel_list_adress, eNR_NBKERNEL_NR);
256
257 for (i = 0; i < eNR_NBKERNEL_NR; i++) {
258 nb_kernel_list_adress[i] = NULL;
259 }
260
261 if (bGenericKernelOnly)
262 {
263 return;
264 }
265
266 nb_kernel_setup_adress(fplog, nb_kernel_list_adress);
267 }
268
269 void do_nonbonded(t_commrec *cr,t_forcerec *fr,
270 rvec x[],rvec f[],t_mdatoms *mdatoms,t_blocka *excl,
271 real egnb[],real egcoul[],real egpol[],rvec box_size,
272 t_nrnb *nrnb,real *lambda, real *dvdl,
273 int nls,int eNL,int flags)
274 {
275 gmx_bool bLR,bDoForces,bForeignLambda;
276 t_nblist * nlist;
277 real * fshift;
278 int n,n0,n1,i,i0,i1,nrnb_ind,sz;
279 t_nblists *nblists;
280 gmx_bool bWater;
281 FILE * fp;
282 int fac=0;
283 int nthreads = 1;
284 int tabletype;
285 int outeriter,inneriter;
286 real * tabledata = NULL;
287 gmx_gbdata_t gbdata;
288
289 nb_kernel_t *kernelptr=NULL;
290 nb_adress_kernel_t *adresskernelptr=NULL;
291
292 gmx_bool bCG; /* for AdresS */
293 int k;/* for AdresS */
294
295 bLR = (flags & GMX_DONB_LR);
296 bDoForces = (flags & GMX_DONB_FORCES);
297 bForeignLambda = (flags & GMX_DONB_FOREIGNLAMBDA);
298
299 bCG = FALSE; /* for AdresS */
300
301 gbdata.gb_epsilon_solvent = fr->gb_epsilon_solvent;
302 gbdata.epsilon_r = fr->epsilon_r;
303 gbdata.gpol = egpol;
304
305 if (!fr->adress_type==eAdressOff && !bDoForces){
306 gmx_fatal(FARGS,"No force kernels not implemeted for adress");
307 }
308
309 if(fr->bAllvsAll)
310 {
311 if(fr->bGB)
312 {
313 #if 0 && defined (GMX_X86_SSE2)
314 # ifdef GMX_DOUBLE
315 if(fr->use_cpu_acceleration)
316 {
317 nb_kernel_allvsallgb_sse2_double(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
318 &outeriter,&inneriter,&fr->AllvsAll_work);
319 }
320 else
321 {
322 nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
323 &outeriter,&inneriter,&fr->AllvsAll_work);
324 }
325 # else /* not double */
326 if(fr->use_cpu_acceleration)
327 {
328 nb_kernel_allvsallgb_sse2_single(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
329 &outeriter,&inneriter,&fr->AllvsAll_work);
330 }
331 else
332 {
333 nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
334 &outeriter,&inneriter,&fr->AllvsAll_work);
335 }
336 # endif /* double/single alt. */
337 #else /* no SSE support compiled in */
338 nb_kernel_allvsallgb(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
339 &outeriter,&inneriter,&fr->AllvsAll_work);
340 #endif
341 inc_nrnb(nrnb,eNR_NBKERNEL_ALLVSALLGB,inneriter);
342 }
343 else
344 {
345 #if 0 && defined (GMX_X86_SSE2)
346 # ifdef GMX_DOUBLE
347 if(fr->use_cpu_acceleration)
348 {
349 nb_kernel_allvsall_sse2_double(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
350 &outeriter,&inneriter,&fr->AllvsAll_work);
351 }
352 else
353 {
354 nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
355 &outeriter,&inneriter,&fr->AllvsAll_work);
356 }
357
358 # else /* not double */
359 if(fr->use_cpu_acceleration)
360 {
361 nb_kernel_allvsall_sse2_single(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
362 &outeriter,&inneriter,&fr->AllvsAll_work);
363 }
364 else
365 {
366 nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
367 &outeriter,&inneriter,&fr->AllvsAll_work);
368 }
369
370 # endif /* double/single check */
371 #else /* No SSE2 support compiled in */
372 nb_kernel_allvsall(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,
373 &outeriter,&inneriter,&fr->AllvsAll_work);
374 #endif
375
376 inc_nrnb(nrnb,eNR_NBKERNEL_ALLVSALL,inneriter);
377 }
378 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,outeriter);
379 return;
380 }
381
382 if (eNL >= 0)
383 {
384 i0 = eNL;
385 i1 = i0+1;
386 }
387 else
388 {
389 i0 = 0;
390 i1 = eNL_NR;
391 }
392
393 if (nls >= 0)
394 {
395 n0 = nls;
396 n1 = nls+1;
397 }
398 else
399 {
400 n0 = 0;
401 n1 = fr->nnblists;
402 }
403
404 if(nb_kernel_list == NULL)
405 {
406 gmx_fatal(FARGS,"gmx_setup_kernels has not been called");
407 }
408
409 fshift = fr->fshift[0];
410
411 for(n=n0; (n<n1); n++)
412 {
413 nblists = &fr->nblists[n];
414 for(i=i0; (i<i1); i++)
415 {
416 outeriter = inneriter = 0;
417
418 if (bLR)
419 {
420 nlist = &(nblists->nlist_lr[i]);
421 }
422 else
423 {
424 nlist = &(nblists->nlist_sr[i]);
425 }
426
427 if (nlist->nri > 0)
428 {
429 nrnb_ind = nlist->il_code;
430
431 if(nrnb_ind==eNR_NBKERNEL_FREE_ENERGY)
432 {
433 /* generic free energy, use combined table */
434 tabledata = nblists->tab.tab;
435 }
436 else
437 {
438 if (bForeignLambda)
439 {
440 /* We don't need the non-perturbed interactions */
441 continue;
442 }
443
444 tabletype = nb_kernel_table[nrnb_ind];
445
446 /* normal kernels, not free energy */
447 if (!bDoForces)
448 {
449 nrnb_ind += eNR_NBKERNEL_NR/2;
450 }
451
452 if(tabletype == TABLE_COMBINED)
453 {
454 tabledata = nblists->tab.tab;
455 }
456 else if(tabletype == TABLE_COUL)
457 {
458 tabledata = nblists->coultab;
459 }
460 else if(tabletype == TABLE_VDW)
461 {
462 tabledata = nblists->vdwtab;
463 }
464 else
465 {
466 tabledata = NULL;
467 }
468 }
469
470 nlist->count = 0;
471
472 if(nlist->free_energy)
473 {
474 if(nlist->ivdw==enbvdwBHAM)
475 {
476 gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
477 }
478
479 gmx_nb_free_energy_kernel(nlist->icoul,
480 nlist->ivdw,
481 nlist->nri,
482 nlist->iinr,
483 nlist->jindex,
484 nlist->jjnr,
485 nlist->shift,
486 fr->shift_vec[0],
487 fshift,
488 nlist->gid,
489 x[0],
490 f[0],
491 mdatoms->chargeA,
492 mdatoms->chargeB,
493 fr->epsfac,
494 fr->k_rf,
495 fr->c_rf,
496 fr->ewaldcoeff,
497 egcoul,
498 mdatoms->typeA,
499 mdatoms->typeB,
500 fr->ntype,
501 fr->nbfp,
502 egnb,
503 nblists->tab.scale,
504 tabledata,
505 lambda[efptCOUL],
506 lambda[efptVDW],
507 dvdl,
508 fr->sc_alphacoul,
509 fr->sc_alphavdw,
510 fr->sc_power,
511 fr->sc_r_power,
512 fr->sc_sigma6_def,
513 fr->sc_sigma6_min,
514 bDoForces,
515 &outeriter,
516 &inneriter);
517 }
518 else if (nlist->enlist == enlistCG_CG)
519 {
520 if (fr->adress_type==eAdressOff){
521 /* Call the charge group based inner loop */
522 gmx_nb_generic_cg_kernel(nlist,
523 fr,
524 mdatoms,
525 x[0],
526 f[0],
527 fshift,
528 egcoul,
529 egnb,
530 nblists->tab.scale,
531 tabledata,
532 &outeriter,
533 &inneriter);
534 }
535 else
536 {
537 /*gmx_nb_generic_adress_kernel(nlist,
538 fr,
539 mdatoms,
540 x[0],
541 f[0],
542 fshift,
543 egcoul,
544 egnb,
545 nblists->tab.scale,
546 tabledata,
547 &outeriter,
548 &inneriter);*/
549 gmx_fatal(FARGS,"Death & horror! Adress cgcg kernel not implemented anymore.\n");
550
551 }
552 }
553 else
554 {
555 /* AdresS*/
556 /* for adress we need to determine for each energy group wether it is explicit or coarse-grained */
557 if (!fr->adress_type == eAdressOff) {
558 bCG = FALSE;
559 if ( !fr->adress_group_explicit[ mdatoms->cENER[nlist->iinr[0]] ] ){
560 bCG=TRUE;
561 }
562 /* If this processor has only explicit atoms (w=1)
563 skip the coarse grained force calculation. Same for
564 only coarsegrained atoms and explicit interactions.
565 Last condition is to make sure that generic kernel is not
566 skipped*/
567 if (mdatoms->pureex && bCG && nb_kernel_list[nrnb_ind] != NULL) continue;
568 if (mdatoms->purecg && !bCG && nb_kernel_list[nrnb_ind] != NULL) continue;
569 }
570
571 if (fr->adress_type == eAdressOff ||
572 mdatoms->pureex ||
573 mdatoms->purecg){
574 /* if we only have to calculate pure cg/ex interactions
575 we can use the faster standard gromacs kernels*/
576 kernelptr = nb_kernel_list[nrnb_ind];
577 }else{
578 /* This processor has hybrid interactions which means
579 * we have to
580 * use our own kernels. We have two kernel types: one that
581 * calculates the forces with the explicit prefactor w1*w2
582 * and one for coarse-grained with (1-w1*w2)
583 * explicit kernels are the second part of the kernel
584 * list */
585 if (!bCG) nrnb_ind += eNR_NBKERNEL_NR/2;
586 adresskernelptr = nb_kernel_list_adress[nrnb_ind];
587 }
588
589 if (kernelptr == NULL && adresskernelptr == NULL)
590 {
591 /* Call a generic nonbonded kernel */
592
593 /* If you want to hack/test your own interactions,
594 * do it in this routine and make sure it is called
595 * by setting the environment variable GMX_NB_GENERIC.
596 */
597 if (fr->adress_type==eAdressOff){
598
599 gmx_nb_generic_kernel(nlist,
600 fr,
601 mdatoms,
602 x[0],
603 f[0],
604 fshift,
605 egcoul,
606 egnb,
607 nblists->tab.scale,
608 tabledata,
609 &outeriter,
610 &inneriter);
611 }else /* do generic AdResS kernels (slow)*/
612 {
613
614 gmx_nb_generic_adress_kernel(nlist,
615 fr,
616 mdatoms,
617 x[0],
618 f[0],
619 fshift,
620 egcoul,
621 egnb,
622 nblists->tab.scale,
623 tabledata,
624 &outeriter,
625 &inneriter,
626 bCG);
627 }
628
629
630 }
631 else
632 {
633 /* Call nonbonded kernel from function pointer */
634 if (kernelptr!=NULL){
635 (*kernelptr)( &(nlist->nri),
636 nlist->iinr,
637 nlist->jindex,
638 nlist->jjnr,
639 nlist->shift,
640 fr->shift_vec[0],
641 fshift,
642 nlist->gid,
643 x[0],
644 f[0],
645 mdatoms->chargeA,
646 &(fr->epsfac),
647 &(fr->k_rf),
648 &(fr->c_rf),
649 egcoul,
650 mdatoms->typeA,
651 &(fr->ntype),
652 fr->nbfp,
653 egnb,
654 &(nblists->tab.scale),
655 tabledata,
656 fr->invsqrta,
657 fr->dvda,
658 &(fr->gbtabscale),
659 fr->gbtab.tab,
660 &nthreads,
661 &(nlist->count),
662 nlist->mtx,
663 &outeriter,
664 &inneriter,
665 (real *)&gbdata);
666 }else if (adresskernelptr != NULL)
667 { /* Adress kernels */
668 (*adresskernelptr)( &(nlist->nri),
669 nlist->iinr,
670 nlist->jindex,
671 nlist->jjnr,
672 nlist->shift,
673 fr->shift_vec[0],
674 fshift,
675 nlist->gid,
676 x[0],
677 f[0],
678 mdatoms->chargeA,
679 &(fr->epsfac),
680 &(fr->k_rf),
681 &(fr->c_rf),
682 egcoul,
683 mdatoms->typeA,
684 &(fr->ntype),
685 fr->nbfp,
686 egnb,
687 &(nblists->tab.scale),
688 tabledata,
689 fr->invsqrta,
690 fr->dvda,
691 &(fr->gbtabscale),
692 fr->gbtab.tab,
693 &nthreads,
694 &(nlist->count),
695 nlist->mtx,
696 &outeriter,
697 &inneriter,
698 fr->adress_ex_forcecap,
699 mdatoms->wf);
700 }
701 }
702 }
703
704 /* Update flop accounting */
705
706 /* Outer loop in kernel */
707 switch (nlist->enlist) {
708 case enlistATOM_ATOM: fac = 1; break;
709 case enlistSPC_ATOM: fac = 3; break;
710 case enlistSPC_SPC: fac = 9; break;
711 case enlistTIP4P_ATOM: fac = 4; break;
712 case enlistTIP4P_TIP4P: fac = 16; break;
713 case enlistCG_CG: fac = 1; break;
714 }
715 inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,fac*outeriter);
716
717 /* inner loop in kernel */
718 inc_nrnb(nrnb,nrnb_ind,inneriter);
719 }
720 }
721 }
722 }
723
724
725 real
726 do_listed_vdw_q(int ftype,int nbonds,
727 const t_iatom iatoms[],const t_iparams iparams[],
728 const rvec x[],rvec f[],rvec fshift[],
729 const t_pbc *pbc,const t_graph *g,
730 real *lambda, real *dvdl,
731 const t_mdatoms *md,
732 const t_forcerec *fr,gmx_grppairener_t *grppener,
733 int *global_atom_index)
734 {
735 static gmx_bool bWarn=FALSE;
736 real eps,r2,*tab,rtab2=0;
737 rvec dx,x14[2],f14[2];
738 int i,ai,aj,itype;
739 int typeA[2]={0,0},typeB[2]={0,1};
740 real chargeA[2]={0,0},chargeB[2];
741 int gid,shift_vir,shift_f;
742 int j_index[] = { 0, 1 };
743 int i0=0,i1=1,i2=2;
744 ivec dt;
745 int outeriter,inneriter;
746 int nthreads = 1;
747 int count;
748 real krf,crf,tabscale;
749 int ntype=0;
750 real *nbfp=NULL;
751 real *egnb=NULL,*egcoul=NULL;
752 t_nblist tmplist;
753 int icoul,ivdw;
754 gmx_bool bMolPBC,bFreeEnergy;
755
756 gmx_bool bCG; /* AdResS*/
757 real wf14[2]={0,0}; /* AdResS*/
758
759 #if GMX_THREAD_SHM_FDECOMP
760 pthread_mutex_t mtx;
761 #else
762 void * mtx = NULL;
763 #endif
764
765
766 #if GMX_THREAD_SHM_FDECOMP
767 pthread_mutex_initialize(&mtx);
768 #endif
769
770 bMolPBC = fr->bMolPBC;
771
772 switch (ftype) {
773 case F_LJ14:
774 case F_LJC14_Q:
775 eps = fr->epsfac*fr->fudgeQQ;
776 ntype = 1;
777 egnb = grppener->ener[egLJ14];
778 egcoul = grppener->ener[egCOUL14];
779 break;
780 case F_LJC_PAIRS_NB:
781 eps = fr->epsfac;
782 ntype = 1;
783 egnb = grppener->ener[egLJSR];
784 egcoul = grppener->ener[egCOULSR];
785 break;
786 default:
787 gmx_fatal(FARGS,"Unknown function type %d in do_nonbonded14",
788 ftype);
789 }
790 tab = fr->tab14.tab;
791 rtab2 = sqr(fr->tab14.r);
792 tabscale = fr->tab14.scale;
793
794 krf = fr->k_rf;
795 crf = fr->c_rf;
796
797 /* Determine the values for icoul/ivdw. */
798 if (fr->bEwald) {
799 icoul = enbcoulOOR;
800 }
801 else if(fr->bcoultab)
802 {
803 icoul = enbcoulTAB;
804 }
805 else if(fr->eeltype == eelRF_NEC)
806 {
807 icoul = enbcoulRF;
808 }
809 else
810 {
811 icoul = enbcoulOOR;
812 }
813
814 if(fr->bvdwtab)
815 {
816 ivdw = enbvdwTAB;
817 }
818 else
819 {
820 ivdw = enbvdwLJ;
821 }
822
823
824 bCG = FALSE; /*Adres*/
825 /* We don't do SSE or altivec here, due to large overhead for 4-fold
826 * unrolling on short lists
827 */
828
829 bFreeEnergy = FALSE;
830 for(i=0; (i<nbonds); )
831 {
832 itype = iatoms[i++];
833 ai = iatoms[i++];
834 aj = iatoms[i++];
835 gid = GID(md->cENER[ai],md->cENER[aj],md->nenergrp);
836
837 if (!fr->adress_type == eAdressOff) {
838 if (fr->adress_group_explicit[md->cENER[ai]] != fr->adress_group_explicit[md->cENER[aj]]){
839 /*exclude cg-ex interaction*/
840 continue;
841 }
842 bCG = !fr->adress_group_explicit[md->cENER[ai]];
843 wf14[0] = md->wf[ai];
844 wf14[1] = md->wf[aj];
845 }
846 switch (ftype) {
847 case F_LJ14:
848 bFreeEnergy =
849 (fr->efep != efepNO &&
850 ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
851 iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
852 iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
853 chargeA[0] = md->chargeA[ai];
854 chargeA[1] = md->chargeA[aj];
855 nbfp = (real *)&(iparams[itype].lj14.c6A);
856 break;
857 case F_LJC14_Q:
858 eps = fr->epsfac*iparams[itype].ljc14.fqq;
859 chargeA[0] = iparams[itype].ljc14.qi;
860 chargeA[1] = iparams[itype].ljc14.qj;
861 nbfp = (real *)&(iparams[itype].ljc14.c6);
862 break;
863 case F_LJC_PAIRS_NB:
864 chargeA[0] = iparams[itype].ljcnb.qi;
865 chargeA[1] = iparams[itype].ljcnb.qj;
866 nbfp = (real *)&(iparams[itype].ljcnb.c6);
867 break;
868 }
869
870 if (!bMolPBC)
871 {
872 /* This is a bonded interaction, atoms are in the same box */
873 shift_f = CENTRAL;
874 r2 = distance2(x[ai],x[aj]);
875 }
876 else
877 {
878 /* Apply full periodic boundary conditions */
879 shift_f = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
880 r2 = norm2(dx);
881 }
882
883 if (r2 >= rtab2)
884 {
885 if (!bWarn)
886 {
887 fprintf(stderr,"Warning: 1-4 interaction between %d and %d "
888 "at distance %.3f which is larger than the 1-4 table size %.3f nm\n",
889 glatnr(global_atom_index,ai),
890 glatnr(global_atom_index,aj),
891 sqrt(r2), sqrt(rtab2));
892 fprintf(stderr,"These are ignored for the rest of the simulation\n");
893 fprintf(stderr,"This usually means your system is exploding,\n"
894 "if not, you should increase table-extension in your mdp file\n"
895 "or with user tables increase the table size\n");
896 bWarn = TRUE;
897 }
898 if (debug)
899 fprintf(debug,"%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
900 x[ai][XX],x[ai][YY],x[ai][ZZ],
901 x[aj][XX],x[aj][YY],x[aj][ZZ],
902 glatnr(global_atom_index,ai),
903 glatnr(global_atom_index,aj),
904 sqrt(r2));
905 }
906 else
907 {
908 copy_rvec(x[ai],x14[0]);
909 copy_rvec(x[aj],x14[1]);
910 clear_rvec(f14[0]);
911 clear_rvec(f14[1]);
912 #ifdef DEBUG
913 fprintf(debug,"LJ14: grp-i=%2d, grp-j=%2d, ngrp=%2d, GID=%d\n",
914 md->cENER[ai],md->cENER[aj],md->nenergrp,gid);
915 #endif
916
917 outeriter = inneriter = count = 0;
918 if (bFreeEnergy)
919 {
920 chargeB[0] = md->chargeB[ai];
921 chargeB[1] = md->chargeB[aj];
922 /* We pass &(iparams[itype].lj14.c6A) as LJ parameter matrix
923 * to the innerloops.
924 * Here we use that the LJ-14 parameters are stored in iparams
925 * as c6A,c12A,c6B,c12B, which are referenced correctly
926 * in the innerloops if we assign type combinations 0-0 and 0-1
927 * to atom pair ai-aj in topologies A and B respectively.
928 */
929
930 /* need to do a bit of a kludge here -- the way it is set up,
931 if the charges change, but the vdw do not, then even though bFreeEnergy is on,
932 it won't work, because all the bonds are perturbed.
933 */
934 if(ivdw==enbvdwBHAM)
935 {
936 gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
937 }
938 count = 0;
939 gmx_nb_free_energy_kernel(icoul,
940 ivdw,
941 i1,
942 &i0,
943 j_index,
944 &i1,
945 &shift_f,
946 fr->shift_vec[0],
947 fshift[0],
948 &gid,
949 x14[0],
950 f14[0],
951 chargeA,
952 chargeB,
953 eps,
954 krf,
955 crf,
956 fr->ewaldcoeff,
957 egcoul,
958 typeA,
959 typeB,
960 ntype,
961 nbfp,
962 egnb,
963 tabscale,
964 tab,
965 lambda[efptCOUL],
966 lambda[efptVDW],
967 dvdl,
968 fr->sc_alphacoul,
969 fr->sc_alphavdw,
970 fr->sc_power,
971 6.0, /* for 1-4's use the 6 power always - 48 power too high because of where they are forced to be */
972 fr->sc_sigma6_def,
973 fr->sc_sigma6_min,
974 TRUE,
975 &outeriter,
976 &inneriter);
977 }
978 else
979 {
980 if (fr->adress_type==eAdressOff || !fr->adress_do_hybridpairs){
981 /* Not perturbed - call kernel 330 */
982 nb_kernel330
983 ( &i1,
984 &i0,
985 j_index,
986 &i1,
987 &shift_f,
988 fr->shift_vec[0],
989 fshift[0],
990 &gid,
991 x14[0],
992 f14[0],
993 chargeA,
994 &eps,
995 &krf,
996 &crf,
997 egcoul,
998 typeA,
999 &ntype,
1000 nbfp,
1001 egnb,
1002 &tabscale,
1003 tab,
1004 NULL,
1005 NULL,
1006 NULL,
1007 NULL,
1008 &nthreads,
1009 &count,
1010 (void *)&mtx,
1011 &outeriter,
1012 &inneriter,
1013 NULL);
1014 } else {
1015 if (bCG) {
1016 nb_kernel330_adress_cg(&i1,
1017 &i0,
1018 j_index,
1019 &i1,
1020 &shift_f,
1021 fr->shift_vec[0],
1022 fshift[0],
1023 &gid,
1024 x14[0],
1025 f14[0],
1026 chargeA,
1027 &eps,
1028 &krf,
1029 &crf,
1030 egcoul,
1031 typeA,
1032 &ntype,
1033 nbfp,
1034 egnb,
1035 &tabscale,
1036 tab,
1037 NULL,
1038 NULL,
1039 NULL,
1040 NULL,
1041 &nthreads,
1042 &count,
1043 (void *) &mtx,
1044 &outeriter,
1045 &inneriter,
1046 fr->adress_ex_forcecap,
1047 wf14);
1048 } else {
1049 nb_kernel330_adress_ex(&i1,
1050 &i0,
1051 j_index,
1052 &i1,
1053 &shift_f,
1054 fr->shift_vec[0],
1055 fshift[0],
1056 &gid,
1057 x14[0],
1058 f14[0],
1059 chargeA,
1060 &eps,
1061 &krf,
1062 &crf,
1063 egcoul,
1064 typeA,
1065 &ntype,
1066 nbfp,
1067 egnb,
1068 &tabscale,
1069 tab,
1070 NULL,
1071 NULL,
1072 NULL,
1073 NULL,
1074 &nthreads,
1075 &count,
1076 (void *) &mtx,
1077 &outeriter,
1078 &inneriter,
1079 fr->adress_ex_forcecap,
1080 wf14);
1081 }
1082
1083 }
1084 }
1085
1086 /* Add the forces */
1087 rvec_inc(f[ai],f14[0]);
1088 rvec_dec(f[aj],f14[0]);
1089
1090 if (g)
1091 {
1092 /* Correct the shift forces using the graph */
1093 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1094 shift_vir = IVEC2IS(dt);
1095 rvec_inc(fshift[shift_vir],f14[0]);
1096 rvec_dec(fshift[CENTRAL],f14[0]);
1097 }
1098
1099 /* flops: eNR_KERNEL_OUTER + eNR_KERNEL330 + 12 */
1100 }
1101 }
1102 return 0.0;
1103 }
1104
1105
The ultimate molecular dynamics simulation package