1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
40 #ifdef GMX_THREAD_SHM_FDECOMP
54 #include "nonbonded.h"
58 #include "gmx_fatal.h"
61 #include "mtop_util.h"
66 /* Uncomment the define below to use the generic charge group - charge group
67 * inner loops (in src/gmxlib/nonbonded/nb_generic_cg.c).
68 * All intra charge-group interaction should be excluded.
69 * There should be no inter charge group exclusions.
70 * There can be no perturbed LJ types or charges.
71 * mdrun currently does NOT check for this.
73 /* #define GMX_CG_INNERLOOP */
77 * E X C L U S I O N H A N D L I N G
81 static void SETEXCL_(t_excl e
[],atom_id i
,atom_id j
)
82 { e
[j
] = e
[j
] | (1<<i
); }
83 static void RMEXCL_(t_excl e
[],atom_id i
,atom_id j
)
84 { e
[j
]=e
[j
] & ~(1<<i
); }
85 static bool ISEXCL_(t_excl e
[],atom_id i
,atom_id j
)
86 { return (bool)(e
[j
] & (1<<i
)); }
87 static bool NOTEXCL_(t_excl e
[],atom_id i
,atom_id j
)
88 { return !(ISEXCL(e
,i
,j
)); }
90 #define SETEXCL(e,i,j) (e)[((atom_id) (j))] |= (1<<((atom_id) (i)))
91 #define RMEXCL(e,i,j) (e)[((atom_id) (j))] &= (~(1<<((atom_id) (i))))
92 #define ISEXCL(e,i,j) (bool) ((e)[((atom_id) (j))] & (1<<((atom_id) (i))))
93 #define NOTEXCL(e,i,j) !(ISEXCL(e,i,j))
96 /************************************************
98 * U T I L I T I E S F O R N S
100 ************************************************/
102 static void reallocate_nblist(t_nblist
*nl
)
106 fprintf(debug
,"reallocating neigborlist il_code=%d, maxnri=%d\n",
107 nl
->il_code
,nl
->maxnri
);
109 srenew(nl
->iinr
, nl
->maxnri
);
110 if (nl
->enlist
== enlistCG_CG
)
112 srenew(nl
->iinr_end
,nl
->maxnri
);
114 srenew(nl
->gid
, nl
->maxnri
);
115 srenew(nl
->shift
, nl
->maxnri
);
116 srenew(nl
->jindex
, nl
->maxnri
+1);
119 /* ivdw/icoul are used to determine the type of interaction, so we
120 * can set an innerloop index here. The obvious choice for this would have
121 * been the vdwtype/coultype values in the forcerecord, but unfortunately
122 * those types are braindead - for instance both Buckingham and normal
123 * Lennard-Jones use the same value (evdwCUT), and a separate boolean variable
124 * to determine which interaction is used. There is further no special value
125 * for 'no interaction'. For backward compatibility with old TPR files we won't
126 * change this in the 3.x series, so when calling this routine you should use:
128 * icoul=0 no coulomb interaction
129 * icoul=1 cutoff standard coulomb
130 * icoul=2 reaction-field coulomb
131 * icoul=3 tabulated coulomb
133 * ivdw=0 no vdw interaction
134 * ivdw=1 standard L-J interaction
136 * ivdw=3 tabulated vdw.
138 * Kind of ugly, but it works.
140 static void init_nblist(t_nblist
*nl_sr
,t_nblist
*nl_lr
,
143 bool bfree
, int enlist
)
175 nl
= (i
== 0) ? nl_sr
: nl_lr
;
176 homenr
= (i
== 0) ? maxsr
: maxlr
;
178 /* Set coul/vdw in neighborlist, and for the normal loops we determine
179 * an index of which one to call.
183 nl
->free_energy
= bfree
;
187 nl
->enlist
= enlistATOM_ATOM
;
188 nl
->il_code
= eNR_NBKERNEL_FREE_ENERGY
;
194 nn
= inloop
[4*icoul
+ ivdw
];
196 /* solvent loops follow directly after the corresponding
197 * ordinary loops, in the order:
199 * SPC, SPC-SPC, TIP4p, TIP4p-TIP4p
203 case enlistATOM_ATOM
:
204 #ifdef GMX_CG_INNERLOOP
205 nl
->enlist
= enlistCG_CG
;
208 case enlistSPC_ATOM
: nn
+= 1; break;
209 case enlistSPC_SPC
: nn
+= 2; break;
210 case enlistTIP4P_ATOM
: nn
+= 3; break;
211 case enlistTIP4P_TIP4P
: nn
+= 4; break;
218 fprintf(debug
,"Initiating neighbourlist type %d for %s interactions,\nwith %d SR, %d LR atoms.\n",
219 nl
->il_code
,ENLISTTYPE(enlist
),maxsr
,maxlr
);
221 /* maxnri is influenced by the number of shifts (maximum is 8)
222 * and the number of energy groups.
223 * If it is not enough, nl memory will be reallocated during the run.
224 * 4 seems to be a reasonable factor, which only causes reallocation
225 * during runs with tiny and many energygroups.
227 nl
->maxnri
= homenr
*4;
236 reallocate_nblist(nl
);
238 #ifdef GMX_THREAD_SHM_FDECOMP
241 pthread_mutex_init(nl
->mtx
,NULL
);
246 void init_neighbor_list(FILE *log
,t_forcerec
*fr
,int homenr
)
248 /* Make maxlr tunable! (does not seem to be a big difference though)
249 * This parameter determines the number of i particles in a long range
250 * neighbourlist. Too few means many function calls, too many means
253 int maxsr
,maxsr_wat
,maxlr
,maxlr_wat
;
254 int icoul
,icoulf
,ivdw
;
256 int enlist_w
,enlist_ww
;
260 /* maxsr = homenr-fr->nWatMol*3; */
265 gmx_fatal(FARGS
,"%s, %d: Negative number of short range atoms.\n"
266 "Call your Gromacs dealer for assistance.",__FILE__
,__LINE__
);
268 /* This is just for initial allocation, so we do not reallocate
269 * all the nlist arrays many times in a row.
270 * The numbers seem very accurate, but they are uncritical.
272 maxsr_wat
= min(fr
->nWatMol
,(homenr
+2)/3);
276 maxlr_wat
= min(maxsr_wat
,maxlr
);
280 maxlr
= maxlr_wat
= 0;
283 /* Determine the values for icoul/ivdw. */
289 else if (fr
->bcoultab
)
293 else if (EEL_RF(fr
->eeltype
))
315 if (fr
->solvent_opt
== esolTIP4P
) {
316 enlist_w
= enlistTIP4P_ATOM
;
317 enlist_ww
= enlistTIP4P_TIP4P
;
319 enlist_w
= enlistSPC_ATOM
;
320 enlist_ww
= enlistSPC_SPC
;
323 for(i
=0; i
<fr
->nnblists
; i
++)
325 nbl
= &(fr
->nblists
[i
]);
326 init_nblist(&nbl
->nlist_sr
[eNL_VDWQQ
],&nbl
->nlist_lr
[eNL_VDWQQ
],
327 maxsr
,maxlr
,ivdw
,icoul
,FALSE
,enlistATOM_ATOM
);
328 init_nblist(&nbl
->nlist_sr
[eNL_VDW
],&nbl
->nlist_lr
[eNL_VDW
],
329 maxsr
,maxlr
,ivdw
,0,FALSE
,enlistATOM_ATOM
);
330 init_nblist(&nbl
->nlist_sr
[eNL_QQ
],&nbl
->nlist_lr
[eNL_QQ
],
331 maxsr
,maxlr
,0,icoul
,FALSE
,enlistATOM_ATOM
);
332 init_nblist(&nbl
->nlist_sr
[eNL_VDWQQ_WATER
],&nbl
->nlist_lr
[eNL_VDWQQ_WATER
],
333 maxsr_wat
,maxlr_wat
,ivdw
,icoul
, FALSE
,enlist_w
);
334 init_nblist(&nbl
->nlist_sr
[eNL_QQ_WATER
],&nbl
->nlist_lr
[eNL_QQ_WATER
],
335 maxsr_wat
,maxlr_wat
,0,icoul
, FALSE
,enlist_w
);
336 init_nblist(&nbl
->nlist_sr
[eNL_VDWQQ_WATERWATER
],&nbl
->nlist_lr
[eNL_VDWQQ_WATERWATER
],
337 maxsr_wat
,maxlr_wat
,ivdw
,icoul
, FALSE
,enlist_ww
);
338 init_nblist(&nbl
->nlist_sr
[eNL_QQ_WATERWATER
],&nbl
->nlist_lr
[eNL_QQ_WATERWATER
],
339 maxsr_wat
,maxlr_wat
,0,icoul
, FALSE
,enlist_ww
);
341 if (fr
->efep
!= efepNO
)
352 init_nblist(&nbl
->nlist_sr
[eNL_VDWQQ_FREE
],&nbl
->nlist_lr
[eNL_VDWQQ_FREE
],
353 maxsr
,maxlr
,ivdw
,icoulf
,TRUE
,enlistATOM_ATOM
);
354 init_nblist(&nbl
->nlist_sr
[eNL_VDW_FREE
],&nbl
->nlist_lr
[eNL_VDW_FREE
],
355 maxsr
,maxlr
,ivdw
,0,TRUE
,enlistATOM_ATOM
);
356 init_nblist(&nbl
->nlist_sr
[eNL_QQ_FREE
],&nbl
->nlist_lr
[eNL_QQ_FREE
],
357 maxsr
,maxlr
,0,icoulf
,TRUE
,enlistATOM_ATOM
);
361 if (fr
->bQMMM
&& fr
->qr
->QMMMscheme
!= eQMMMschemeoniom
)
363 init_nblist(&fr
->QMMMlist_sr
,&fr
->QMMMlist_lr
,
364 maxsr
,maxlr
,0,icoul
,FALSE
,enlistATOM_ATOM
);
367 fr
->ns
.nblist_initialized
=TRUE
;
370 static void reset_nblist(t_nblist
*nl
)
381 static void reset_neighbor_list(t_forcerec
*fr
,bool bLR
,int nls
,int eNL
)
387 reset_nblist(&(fr
->nblists
[nls
].nlist_lr
[eNL
]));
391 for(n
=0; n
<fr
->nnblists
; n
++)
393 for(i
=0; i
<eNL_NR
; i
++)
395 reset_nblist(&(fr
->nblists
[n
].nlist_sr
[i
]));
400 /* only reset the short-range nblist */
401 reset_nblist(&(fr
->QMMMlist_sr
));
409 static inline void new_i_nblist(t_nblist
*nlist
,
410 bool bLR
,atom_id i_atom
,int shift
,int gid
)
416 /* Check whether we have to increase the i counter */
418 (nlist
->iinr
[nri
] != i_atom
) ||
419 (nlist
->shift
[nri
] != shift
) ||
420 (nlist
->gid
[nri
] != gid
))
422 /* This is something else. Now see if any entries have
423 * been added in the list of the previous atom.
426 ((nlist
->jindex
[nri
+1] > nlist
->jindex
[nri
]) &&
427 (nlist
->gid
[nri
] != -1)))
429 /* If so increase the counter */
432 if (nlist
->nri
>= nlist
->maxnri
)
434 nlist
->maxnri
+= over_alloc_large(nlist
->nri
);
435 reallocate_nblist(nlist
);
438 /* Set the number of neighbours and the atom number */
439 nlist
->jindex
[nri
+1] = nlist
->jindex
[nri
];
440 nlist
->iinr
[nri
] = i_atom
;
441 nlist
->gid
[nri
] = gid
;
442 nlist
->shift
[nri
] = shift
;
446 static inline void close_i_nblist(t_nblist
*nlist
)
448 int nri
= nlist
->nri
;
453 nlist
->jindex
[nri
+1] = nlist
->nrj
;
455 len
=nlist
->nrj
- nlist
->jindex
[nri
];
457 /* nlist length for water i molecules is treated statically
460 if (len
> nlist
->maxlen
)
467 static inline void close_nblist(t_nblist
*nlist
)
469 /* Only close this nblist when it has been initialized */
476 static inline void close_neighbor_list(t_forcerec
*fr
,bool bLR
,int nls
,int eNL
,
477 bool bMakeQMMMnblist
)
481 if (bMakeQMMMnblist
) {
484 close_nblist(&(fr
->QMMMlist_lr
));
488 close_nblist(&(fr
->QMMMlist_sr
));
495 close_nblist(&(fr
->nblists
[nls
].nlist_lr
[eNL
]));
499 for(n
=0; n
<fr
->nnblists
; n
++)
501 for(i
=0; (i
<eNL_NR
); i
++)
503 close_nblist(&(fr
->nblists
[n
].nlist_sr
[i
]));
510 static inline void add_j_to_nblist(t_nblist
*nlist
,atom_id j_atom
,bool bLR
)
514 if (nlist
->nrj
>= nlist
->maxnrj
)
516 nlist
->maxnrj
= over_alloc_small(nlist
->nrj
+ 1);
518 fprintf(debug
,"Increasing %s nblist %s j size to %d\n",
519 bLR
? "LR" : "SR",nrnb_str(nlist
->il_code
),nlist
->maxnrj
);
521 srenew(nlist
->jjnr
,nlist
->maxnrj
);
524 nlist
->jjnr
[nrj
] = j_atom
;
528 static inline void add_j_to_nblist_cg(t_nblist
*nlist
,
529 atom_id j_start
,int j_end
,bool bLR
)
533 if (nlist
->nrj
>= nlist
->maxnrj
)
535 nlist
->maxnrj
= over_alloc_small(nlist
->nrj
+ 1);
537 fprintf(debug
,"Increasing %s nblist %s j size to %d\n",
538 bLR
? "LR" : "SR",nrnb_str(nlist
->il_code
),nlist
->maxnrj
);
540 srenew(nlist
->jjnr
,nlist
->maxnrj
);
541 srenew(nlist
->jjnr_end
,nlist
->maxnrj
);
544 nlist
->jjnr
[nrj
] = j_start
;
545 nlist
->jjnr_end
[nrj
] = j_end
;
549 #ifndef GMX_CG_INNERLOOP
551 put_in_list(bool bHaveVdW
[],
567 /* The a[] index has been removed,
568 * to put it back in i_atom should be a[i0] and jj should be a[jj].
573 t_nblist
* vdwc_free
= NULL
;
574 t_nblist
* vdw_free
= NULL
;
575 t_nblist
* coul_free
= NULL
;
576 t_nblist
* vdwc_ww
= NULL
;
577 t_nblist
* coul_ww
= NULL
;
579 int i
,j
,jcg
,igid
,gid
,nbl_ind
,ind_ij
;
580 atom_id jj
,jj0
,jj1
,i_atom
;
585 real
*charge
,*chargeB
;
587 bool bFreeEnergy
,bFree
,bFreeJ
,bNotEx
,*bPert
;
588 bool bDoVdW_i
,bDoCoul_i
,bDoCoul_i_sol
;
592 /* Copy some pointers */
594 charge
= md
->chargeA
;
595 chargeB
= md
->chargeB
;
598 bPert
= md
->bPerturbed
;
602 nicg
= index
[icg
+1]-i0
;
604 /* Get the i charge group info */
605 igid
= GET_CGINFO_GID(cginfo
[icg
]);
606 iwater
= GET_CGINFO_SOLOPT(cginfo
[icg
]);
611 /* Check if any of the particles involved are perturbed.
612 * If not we can do the cheaper normal put_in_list
613 * and use more solvent optimization.
615 for(i
=0; i
<nicg
; i
++)
617 bFreeEnergy
|= bPert
[i0
+i
];
619 /* Loop over the j charge groups */
620 for(j
=0; (j
<nj
&& !bFreeEnergy
); j
++)
625 /* Finally loop over the atoms in the j-charge group */
626 for(jj
=jj0
; jj
<jj1
; jj
++)
628 bFreeEnergy
|= bPert
[jj
];
633 /* Unpack pointers to neighbourlist structs */
634 if (fr
->nnblists
== 1)
640 nbl_ind
= fr
->gid2nblists
[GID(igid
,jgid
,ngid
)];
644 nlist
= fr
->nblists
[nbl_ind
].nlist_lr
;
648 nlist
= fr
->nblists
[nbl_ind
].nlist_sr
;
651 if (iwater
!= esolNO
)
653 vdwc
= &nlist
[eNL_VDWQQ_WATER
];
654 vdw
= &nlist
[eNL_VDW
];
655 coul
= &nlist
[eNL_QQ_WATER
];
656 #ifndef DISABLE_WATERWATER_NLIST
657 vdwc_ww
= &nlist
[eNL_VDWQQ_WATERWATER
];
658 coul_ww
= &nlist
[eNL_QQ_WATERWATER
];
665 coul
= (bLR
) ? &fr
->QMMMlist_lr
: &fr
->QMMMlist_sr
;
669 vdwc
= &nlist
[eNL_VDWQQ
];
670 vdw
= &nlist
[eNL_VDW
];
671 coul
= &nlist
[eNL_QQ
];
676 if (iwater
!= esolNO
)
678 /* Loop over the atoms in the i charge group */
680 gid
= GID(igid
,jgid
,ngid
);
681 /* Create new i_atom for each energy group */
682 if (bDoCoul
&& bDoVdW
)
684 new_i_nblist(vdwc
,bLR
,i_atom
,shift
,gid
);
685 #ifndef DISABLE_WATERWATER_NLIST
686 new_i_nblist(vdwc_ww
,bLR
,i_atom
,shift
,gid
);
691 new_i_nblist(vdw
,bLR
,i_atom
,shift
,gid
);
695 new_i_nblist(coul
,bLR
,i_atom
,shift
,gid
);
696 #ifndef DISABLE_WATERWATER_NLIST
697 new_i_nblist(coul_ww
,bLR
,i_atom
,shift
,gid
);
700 /* Loop over the j charge groups */
701 for(j
=0; (j
<nj
); j
++)
711 jwater
= GET_CGINFO_SOLOPT(cginfo
[jcg
]);
713 if (iwater
== esolSPC
&& jwater
== esolSPC
)
715 /* Interaction between two SPC molecules */
718 /* VdW only - only first atoms in each water interact */
719 add_j_to_nblist(vdw
,jj0
,bLR
);
723 #ifdef DISABLE_WATERWATER_NLIST
724 /* Add entries for the three atoms - only do VdW if we need to */
727 add_j_to_nblist(coul
,jj0
,bLR
);
731 add_j_to_nblist(vdwc
,jj0
,bLR
);
733 add_j_to_nblist(coul
,jj0
+1,bLR
);
734 add_j_to_nblist(coul
,jj0
+2,bLR
);
736 /* One entry for the entire water-water interaction */
739 add_j_to_nblist(coul_ww
,jj0
,bLR
);
743 add_j_to_nblist(vdwc_ww
,jj0
,bLR
);
748 else if (iwater
== esolTIP4P
&& jwater
== esolTIP4P
)
750 /* Interaction between two TIP4p molecules */
753 /* VdW only - only first atoms in each water interact */
754 add_j_to_nblist(vdw
,jj0
,bLR
);
758 #ifdef DISABLE_WATERWATER_NLIST
759 /* Add entries for the four atoms - only do VdW if we need to */
762 add_j_to_nblist(vdw
,jj0
,bLR
);
764 add_j_to_nblist(coul
,jj0
+1,bLR
);
765 add_j_to_nblist(coul
,jj0
+2,bLR
);
766 add_j_to_nblist(coul
,jj0
+3,bLR
);
768 /* One entry for the entire water-water interaction */
771 add_j_to_nblist(coul_ww
,jj0
,bLR
);
775 add_j_to_nblist(vdwc_ww
,jj0
,bLR
);
782 /* j charge group is not water, but i is.
783 * Add entries to the water-other_atom lists; the geometry of the water
784 * molecule doesn't matter - that is taken care of in the nonbonded kernel,
785 * so we don't care if it is SPC or TIP4P...
792 for(jj
=jj0
; (jj
<jj1
); jj
++)
796 add_j_to_nblist(coul
,jj
,bLR
);
802 for(jj
=jj0
; (jj
<jj1
); jj
++)
804 if (bHaveVdW
[type
[jj
]])
806 add_j_to_nblist(vdw
,jj
,bLR
);
812 /* _charge_ _groups_ interact with both coulomb and LJ */
813 /* Check which atoms we should add to the lists! */
814 for(jj
=jj0
; (jj
<jj1
); jj
++)
816 if (bHaveVdW
[type
[jj
]])
820 add_j_to_nblist(vdwc
,jj
,bLR
);
824 add_j_to_nblist(vdw
,jj
,bLR
);
827 else if (charge
[jj
] != 0)
829 add_j_to_nblist(coul
,jj
,bLR
);
836 close_i_nblist(coul
);
837 close_i_nblist(vdwc
);
838 #ifndef DISABLE_WATERWATER_NLIST
839 close_i_nblist(coul_ww
);
840 close_i_nblist(vdwc_ww
);
845 /* QMMM atoms as i charge groups */
848 /* Loop over atoms in the ith charge group */
852 gid
= GID(igid
,jgid
,ngid
);
853 /* Create new i_atom for each energy group */
854 new_i_nblist(coul
,bLR
,i_atom
,shift
,gid
);
856 /* Loop over the j charge groups */
861 /* Charge groups cannot have QM and MM atoms simultaneously */
866 /* Finally loop over the atoms in the j-charge group */
867 for(jj
=jj0
; jj
<jj1
; jj
++)
869 bNotEx
= NOTEXCL(bExcl
,i
,jj
);
871 add_j_to_nblist(coul
,jj
,bLR
);
875 close_i_nblist(coul
);
881 /* no solvent as i charge group */
882 /* Loop over the atoms in the i charge group */
883 for(i
=0; i
<nicg
; i
++)
886 gid
= GID(igid
,jgid
,ngid
);
889 /* Create new i_atom for each energy group */
890 if (bDoVdW
&& bDoCoul
)
892 new_i_nblist(vdwc
,bLR
,i_atom
,shift
,gid
);
896 new_i_nblist(vdw
,bLR
,i_atom
,shift
,gid
);
900 new_i_nblist(coul
,bLR
,i_atom
,shift
,gid
);
902 bDoVdW_i
= (bDoVdW
&& bHaveVdW
[type
[i_atom
]]);
903 bDoCoul_i
= (bDoCoul
&& qi
!=0);
905 if (bDoVdW_i
|| bDoCoul_i
)
907 /* Loop over the j charge groups */
908 for(j
=0; (j
<nj
); j
++)
912 /* Check for large charge groups */
923 /* Finally loop over the atoms in the j-charge group */
924 for(jj
=jj0
; jj
<jj1
; jj
++)
926 bNotEx
= NOTEXCL(bExcl
,i
,jj
);
934 add_j_to_nblist(coul
,jj
,bLR
);
939 if (bHaveVdW
[type
[jj
]])
941 add_j_to_nblist(vdw
,jj
,bLR
);
946 if (bHaveVdW
[type
[jj
]])
950 add_j_to_nblist(vdwc
,jj
,bLR
);
954 add_j_to_nblist(vdw
,jj
,bLR
);
957 else if (charge
[jj
] != 0)
959 add_j_to_nblist(coul
,jj
,bLR
);
967 close_i_nblist(coul
);
968 close_i_nblist(vdwc
);
974 /* we are doing free energy */
975 vdwc_free
= &nlist
[eNL_VDWQQ_FREE
];
976 vdw_free
= &nlist
[eNL_VDW_FREE
];
977 coul_free
= &nlist
[eNL_QQ_FREE
];
978 /* Loop over the atoms in the i charge group */
979 for(i
=0; i
<nicg
; i
++)
982 gid
= GID(igid
,jgid
,ngid
);
984 qiB
= chargeB
[i_atom
];
986 /* Create new i_atom for each energy group */
987 if (bDoVdW
&& bDoCoul
)
988 new_i_nblist(vdwc
,bLR
,i_atom
,shift
,gid
);
990 new_i_nblist(vdw
,bLR
,i_atom
,shift
,gid
);
992 new_i_nblist(coul
,bLR
,i_atom
,shift
,gid
);
994 new_i_nblist(vdw_free
,bLR
,i_atom
,shift
,gid
);
995 new_i_nblist(coul_free
,bLR
,i_atom
,shift
,gid
);
996 new_i_nblist(vdwc_free
,bLR
,i_atom
,shift
,gid
);
998 bDoVdW_i
= (bDoVdW
&&
999 (bHaveVdW
[type
[i_atom
]] || bHaveVdW
[typeB
[i_atom
]]));
1000 bDoCoul_i
= (bDoCoul
&& (qi
!=0 || qiB
!=0));
1001 /* For TIP4P the first atom does not have a charge,
1002 * but the last three do. So we should still put an atom
1003 * without LJ but with charge in the water-atom neighborlist
1004 * for a TIP4p i charge group.
1005 * For SPC type water the first atom has LJ and charge,
1006 * so there is no such problem.
1008 if (iwater
== esolNO
)
1010 bDoCoul_i_sol
= bDoCoul_i
;
1014 bDoCoul_i_sol
= bDoCoul
;
1017 if (bDoVdW_i
|| bDoCoul_i_sol
)
1019 /* Loop over the j charge groups */
1020 for(j
=0; (j
<nj
); j
++)
1024 /* Check for large charge groups */
1035 /* Finally loop over the atoms in the j-charge group */
1036 bFree
= bPert
[i_atom
];
1037 for(jj
=jj0
; (jj
<jj1
); jj
++)
1039 bFreeJ
= bFree
|| bPert
[jj
];
1040 /* Complicated if, because the water H's should also
1041 * see perturbed j-particles
1043 if (iwater
==esolNO
|| i
==0 || bFreeJ
)
1045 bNotEx
= NOTEXCL(bExcl
,i
,jj
);
1053 if (charge
[jj
]!=0 || chargeB
[jj
]!=0)
1055 add_j_to_nblist(coul_free
,jj
,bLR
);
1058 else if (!bDoCoul_i
)
1060 if (bHaveVdW
[type
[jj
]] || bHaveVdW
[typeB
[jj
]])
1062 add_j_to_nblist(vdw_free
,jj
,bLR
);
1067 if (bHaveVdW
[type
[jj
]] || bHaveVdW
[typeB
[jj
]])
1069 if (charge
[jj
]!=0 || chargeB
[jj
]!=0)
1071 add_j_to_nblist(vdwc_free
,jj
,bLR
);
1075 add_j_to_nblist(vdw_free
,jj
,bLR
);
1078 else if (charge
[jj
]!=0 || chargeB
[jj
]!=0)
1079 add_j_to_nblist(coul_free
,jj
,bLR
);
1084 /* This is done whether or not bWater is set */
1085 if (charge
[jj
] != 0)
1087 add_j_to_nblist(coul
,jj
,bLR
);
1090 else if (!bDoCoul_i_sol
)
1092 if (bHaveVdW
[type
[jj
]])
1094 add_j_to_nblist(vdw
,jj
,bLR
);
1099 if (bHaveVdW
[type
[jj
]])
1101 if (charge
[jj
] != 0)
1103 add_j_to_nblist(vdwc
,jj
,bLR
);
1107 add_j_to_nblist(vdw
,jj
,bLR
);
1110 else if (charge
[jj
] != 0)
1112 add_j_to_nblist(coul
,jj
,bLR
);
1120 close_i_nblist(vdw
);
1121 close_i_nblist(coul
);
1122 close_i_nblist(vdwc
);
1123 close_i_nblist(vdw_free
);
1124 close_i_nblist(coul_free
);
1125 close_i_nblist(vdwc_free
);
1131 #ifdef GMX_CG_INNERLOOP
1133 put_in_list(bool bHaveVdW
[],
1149 int igid
,gid
,nbl_ind
;
1153 igid
= GET_CGINFO_GID(fr
->cginfo
[icg
]);
1154 gid
= GID(igid
,jgid
,ngid
);
1156 /* Unpack pointers to neighbourlist structs */
1157 if (fr
->nnblists
== 1)
1163 nbl_ind
= fr
->gid2nblists
[gid
];
1167 vdwc
= &fr
->nblists
[nbl_ind
].nlist_lr
[eNL_VDWQQ
];
1171 vdwc
= &fr
->nblists
[nbl_ind
].nlist_sr
[eNL_VDWQQ
];
1174 /* Make a new neighbor list for charge group icg.
1175 * Currently simply one neighbor list is made with LJ and Coulomb.
1176 * If required, zerp interactions could be removed here
1177 * or in the force loop.
1179 new_i_nblist(vdwc
,bLR
,index
[icg
],shift
,gid
);
1180 vdwc
->iinr_end
[vdwc
->nri
] = index
[icg
+1];
1182 for(j
=0; (j
<nj
); j
++)
1187 /* Here we simply add the j charge group jcg to the list
1188 * without checking exclusions, LJ interactions or charges.
1190 add_j_to_nblist_cg(vdwc
,index
[jcg
],index
[jcg
+1],bLR
);
1194 close_i_nblist(vdwc
);
1198 static void setexcl(atom_id start
,atom_id end
,t_blocka
*excl
,bool b
,
1205 for(i
=start
; i
<end
; i
++)
1207 for(k
=excl
->index
[i
]; k
<excl
->index
[i
+1]; k
++)
1209 SETEXCL(bexcl
,i
-start
,excl
->a
[k
]);
1215 for(i
=start
; i
<end
; i
++)
1217 for(k
=excl
->index
[i
]; k
<excl
->index
[i
+1]; k
++)
1219 RMEXCL(bexcl
,i
-start
,excl
->a
[k
]);
1225 int calc_naaj(int icg
,int cgtot
)
1229 if ((cgtot
% 2) == 1)
1231 /* Odd number of charge groups, easy */
1232 naaj
= 1 + (cgtot
/2);
1234 else if ((cgtot
% 4) == 0)
1236 /* Multiple of four is hard */
1273 fprintf(log
,"naaj=%d\n",naaj
);
1279 /************************************************
1281 * S I M P L E C O R E S T U F F
1283 ************************************************/
1285 static real
calc_image_tric(rvec xi
,rvec xj
,matrix box
,
1286 rvec b_inv
,int *shift
)
1288 /* This code assumes that the cut-off is smaller than
1289 * a half times the smallest diagonal element of the box.
1296 /* Compute diff vector */
1297 dz
= xj
[ZZ
] - xi
[ZZ
];
1298 dy
= xj
[YY
] - xi
[YY
];
1299 dx
= xj
[XX
] - xi
[XX
];
1301 /* Perform NINT operation, using trunc operation, therefore
1302 * we first add 2.5 then subtract 2 again
1304 tz
= dz
*b_inv
[ZZ
] + h25
;
1306 dz
-= tz
*box
[ZZ
][ZZ
];
1307 dy
-= tz
*box
[ZZ
][YY
];
1308 dx
-= tz
*box
[ZZ
][XX
];
1310 ty
= dy
*b_inv
[YY
] + h25
;
1312 dy
-= ty
*box
[YY
][YY
];
1313 dx
-= ty
*box
[YY
][XX
];
1315 tx
= dx
*b_inv
[XX
]+h25
;
1317 dx
-= tx
*box
[XX
][XX
];
1319 /* Distance squared */
1320 r2
= (dx
*dx
) + (dy
*dy
) + (dz
*dz
);
1322 *shift
= XYZ2IS(tx
,ty
,tz
);
1327 static real
calc_image_rect(rvec xi
,rvec xj
,rvec box_size
,
1328 rvec b_inv
,int *shift
)
1336 /* Compute diff vector */
1337 dx
= xj
[XX
] - xi
[XX
];
1338 dy
= xj
[YY
] - xi
[YY
];
1339 dz
= xj
[ZZ
] - xi
[ZZ
];
1341 /* Perform NINT operation, using trunc operation, therefore
1342 * we first add 1.5 then subtract 1 again
1344 tx
= dx
*b_inv
[XX
] + h15
;
1345 ty
= dy
*b_inv
[YY
] + h15
;
1346 tz
= dz
*b_inv
[ZZ
] + h15
;
1351 /* Correct diff vector for translation */
1352 ddx
= tx
*box_size
[XX
] - dx
;
1353 ddy
= ty
*box_size
[YY
] - dy
;
1354 ddz
= tz
*box_size
[ZZ
] - dz
;
1356 /* Distance squared */
1357 r2
= (ddx
*ddx
) + (ddy
*ddy
) + (ddz
*ddz
);
1359 *shift
= XYZ2IS(tx
,ty
,tz
);
1364 static void add_simple(t_ns_buf
*nsbuf
,int nrj
,atom_id cg_j
,
1365 bool bHaveVdW
[],int ngid
,t_mdatoms
*md
,
1366 int icg
,int jgid
,t_block
*cgs
,t_excl bexcl
[],
1367 int shift
,t_forcerec
*fr
)
1369 if (nsbuf
->nj
+ nrj
> MAX_CG
)
1371 put_in_list(bHaveVdW
,ngid
,md
,icg
,jgid
,nsbuf
->ncg
,nsbuf
->jcg
,
1372 cgs
->index
,bexcl
,shift
,fr
,FALSE
,TRUE
,TRUE
,FALSE
);
1373 /* Reset buffer contents */
1374 nsbuf
->ncg
= nsbuf
->nj
= 0;
1376 nsbuf
->jcg
[nsbuf
->ncg
++] = cg_j
;
1380 static void ns_inner_tric(rvec x
[],int icg
,int *i_egp_flags
,
1381 int njcg
,atom_id jcg
[],
1382 matrix box
,rvec b_inv
,real rcut2
,
1383 t_block
*cgs
,t_ns_buf
**ns_buf
,
1384 bool bHaveVdW
[],int ngid
,t_mdatoms
*md
,
1385 t_excl bexcl
[],t_forcerec
*fr
)
1389 int *cginfo
=fr
->cginfo
;
1390 atom_id cg_j
,*cgindex
;
1393 cgindex
= cgs
->index
;
1395 for(j
=0; (j
<njcg
); j
++)
1398 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1399 if (calc_image_tric(x
[icg
],x
[cg_j
],box
,b_inv
,&shift
) < rcut2
)
1401 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1402 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1404 add_simple(&ns_buf
[jgid
][shift
],nrj
,cg_j
,
1405 bHaveVdW
,ngid
,md
,icg
,jgid
,cgs
,bexcl
,shift
,fr
);
1411 static void ns_inner_rect(rvec x
[],int icg
,int *i_egp_flags
,
1412 int njcg
,atom_id jcg
[],
1413 bool bBox
,rvec box_size
,rvec b_inv
,real rcut2
,
1414 t_block
*cgs
,t_ns_buf
**ns_buf
,
1415 bool bHaveVdW
[],int ngid
,t_mdatoms
*md
,
1416 t_excl bexcl
[],t_forcerec
*fr
)
1420 int *cginfo
=fr
->cginfo
;
1421 atom_id cg_j
,*cgindex
;
1424 cgindex
= cgs
->index
;
1428 for(j
=0; (j
<njcg
); j
++)
1431 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1432 if (calc_image_rect(x
[icg
],x
[cg_j
],box_size
,b_inv
,&shift
) < rcut2
)
1434 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1435 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1437 add_simple(&ns_buf
[jgid
][shift
],nrj
,cg_j
,
1438 bHaveVdW
,ngid
,md
,icg
,jgid
,cgs
,bexcl
,shift
,fr
);
1445 for(j
=0; (j
<njcg
); j
++)
1448 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1449 if ((rcut2
== 0) || (distance2(x
[icg
],x
[cg_j
]) < rcut2
)) {
1450 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1451 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1453 add_simple(&ns_buf
[jgid
][CENTRAL
],nrj
,cg_j
,
1454 bHaveVdW
,ngid
,md
,icg
,jgid
,cgs
,bexcl
,CENTRAL
,fr
);
1461 /* ns_simple_core needs to be adapted for QMMM still 2005 */
1463 static int ns_simple_core(t_forcerec
*fr
,
1464 gmx_localtop_t
*top
,
1466 matrix box
,rvec box_size
,
1467 t_excl bexcl
[],atom_id
*aaj
,
1468 int ngid
,t_ns_buf
**ns_buf
,
1473 int nsearch
,icg
,jcg
,igid
,i0
,nri
,nn
;
1476 /* atom_id *i_atoms; */
1477 t_block
*cgs
=&(top
->cgs
);
1478 t_blocka
*excl
=&(top
->excls
);
1481 bool bBox
,bTriclinic
;
1484 rlist2
= sqr(fr
->rlist
);
1486 bBox
= (fr
->ePBC
!= epbcNONE
);
1489 for(m
=0; (m
<DIM
); m
++)
1491 b_inv
[m
] = divide(1.0,box_size
[m
]);
1493 bTriclinic
= TRICLINIC(box
);
1500 cginfo
= fr
->cginfo
;
1503 for (icg
=fr
->cg0
; (icg
<fr
->hcg
); icg
++)
1506 i0 = cgs->index[icg];
1507 nri = cgs->index[icg+1]-i0;
1508 i_atoms = &(cgs->a[i0]);
1509 i_eg_excl = fr->eg_excl + ngid*md->cENER[*i_atoms];
1510 setexcl(nri,i_atoms,excl,TRUE,bexcl);
1512 igid
= GET_CGINFO_GID(cginfo
[icg
]);
1513 i_egp_flags
= fr
->egp_flags
+ ngid
*igid
;
1514 setexcl(cgs
->index
[icg
],cgs
->index
[icg
+1],excl
,TRUE
,bexcl
);
1516 naaj
=calc_naaj(icg
,cgs
->nr
);
1519 ns_inner_tric(fr
->cg_cm
,icg
,i_egp_flags
,naaj
,&(aaj
[icg
]),
1520 box
,b_inv
,rlist2
,cgs
,ns_buf
,
1521 bHaveVdW
,ngid
,md
,bexcl
,fr
);
1525 ns_inner_rect(fr
->cg_cm
,icg
,i_egp_flags
,naaj
,&(aaj
[icg
]),
1526 bBox
,box_size
,b_inv
,rlist2
,cgs
,ns_buf
,
1527 bHaveVdW
,ngid
,md
,bexcl
,fr
);
1531 for(nn
=0; (nn
<ngid
); nn
++)
1533 for(k
=0; (k
<SHIFTS
); k
++)
1535 nsbuf
= &(ns_buf
[nn
][k
]);
1538 put_in_list(bHaveVdW
,ngid
,md
,icg
,nn
,nsbuf
->ncg
,nsbuf
->jcg
,
1539 cgs
->index
,bexcl
,k
,fr
,FALSE
,TRUE
,TRUE
,FALSE
);
1540 nsbuf
->ncg
=nsbuf
->nj
=0;
1544 /* setexcl(nri,i_atoms,excl,FALSE,bexcl); */
1545 setexcl(cgs
->index
[icg
],cgs
->index
[icg
+1],excl
,FALSE
,bexcl
);
1547 close_neighbor_list(fr
,FALSE
,-1,-1,FALSE
);
1552 /************************************************
1554 * N S 5 G R I D S T U F F
1556 ************************************************/
1558 static inline void get_dx(int Nx
,real gridx
,real rc2
,int xgi
,real x
,
1559 int *dx0
,int *dx1
,real
*dcx2
)
1587 for(i
=xgi0
; i
>=0; i
--)
1589 dcx
= (i
+1)*gridx
-x
;
1596 for(i
=xgi1
; i
<Nx
; i
++)
1609 static inline void get_dx_dd(int Nx
,real gridx
,real rc2
,int xgi
,real x
,
1610 int ncpddc
,int shift_min
,int shift_max
,
1611 int *g0
,int *g1
,real
*dcx2
)
1614 int g_min
,g_max
,shift_home
;
1647 g_min
= (shift_min
== shift_home
? 0 : ncpddc
);
1648 g_max
= (shift_max
== shift_home
? ncpddc
- 1 : Nx
- 1);
1655 else if (shift_max
< 0)
1670 /* Check one grid cell down */
1671 dcx
= ((*g0
- 1) + 1)*gridx
- x
;
1683 /* Check one grid cell up */
1684 dcx
= (*g1
+ 1)*gridx
- x
;
1696 #define sqr(x) ((x)*(x))
1697 #define calc_dx2(XI,YI,ZI,y) (sqr(XI-y[XX]) + sqr(YI-y[YY]) + sqr(ZI-y[ZZ]))
1698 #define calc_cyl_dx2(XI,YI,y) (sqr(XI-y[XX]) + sqr(YI-y[YY]))
1699 /****************************************************
1701 * F A S T N E I G H B O R S E A R C H I N G
1703 * Optimized neighboursearching routine using grid
1704 * at least 1x1x1, see GROMACS manual
1706 ****************************************************/
1708 static void do_longrange(t_commrec
*cr
,gmx_localtop_t
*top
,t_forcerec
*fr
,
1709 int ngid
,t_mdatoms
*md
,int icg
,
1711 atom_id lr
[],t_excl bexcl
[],int shift
,
1712 rvec x
[],rvec box_size
,t_nrnb
*nrnb
,
1713 real lambda
,real
*dvdlambda
,
1714 gmx_grppairener_t
*grppener
,
1715 bool bDoVdW
,bool bDoCoul
,
1716 bool bEvaluateNow
,bool bHaveVdW
[],
1717 bool bDoForces
,rvec
*f
)
1722 for(n
=0; n
<fr
->nnblists
; n
++)
1724 for(i
=0; (i
<eNL_NR
); i
++)
1726 nl
= &fr
->nblists
[n
].nlist_lr
[i
];
1727 if ((nl
->nri
> nl
->maxnri
-32) || bEvaluateNow
)
1729 close_neighbor_list(fr
,TRUE
,n
,i
,FALSE
);
1730 /* Evaluate the energies and forces */
1731 do_nonbonded(cr
,fr
,x
,f
,md
,NULL
,
1732 grppener
->ener
[fr
->bBHAM
? egBHAMLR
: egLJLR
],
1733 grppener
->ener
[egCOULLR
],
1734 grppener
->ener
[egGB
],box_size
,
1735 nrnb
,lambda
,dvdlambda
,n
,i
,
1736 GMX_DONB_LR
| GMX_DONB_FORCES
);
1738 reset_neighbor_list(fr
,TRUE
,n
,i
);
1745 /* Put the long range particles in a list */
1746 /* Since do_longrange is never called for QMMM, bQMMM is FALSE */
1747 put_in_list(bHaveVdW
,ngid
,md
,icg
,jgid
,nlr
,lr
,top
->cgs
.index
,
1748 bexcl
,shift
,fr
,TRUE
,bDoVdW
,bDoCoul
,FALSE
);
1752 static int nsgrid_core(FILE *log
,t_commrec
*cr
,t_forcerec
*fr
,
1753 matrix box
,rvec box_size
,int ngid
,
1754 gmx_localtop_t
*top
,
1755 t_grid
*grid
,rvec x
[],
1756 t_excl bexcl
[],bool *bExcludeAlleg
,
1757 t_nrnb
*nrnb
,t_mdatoms
*md
,
1758 real lambda
,real
*dvdlambda
,
1759 gmx_grppairener_t
*grppener
,
1761 bool bDoLongRange
,bool bDoForces
,rvec
*f
,
1762 bool bMakeQMMMnblist
)
1765 atom_id
**nl_lr_ljc
,**nl_lr_one
,**nl_sr
;
1766 int *nlr_ljc
,*nlr_one
,*nsr
;
1767 gmx_domdec_t
*dd
=NULL
;
1768 t_block
*cgs
=&(top
->cgs
);
1769 int *cginfo
=fr
->cginfo
;
1770 /* atom_id *i_atoms,*cgsindex=cgs->index; */
1772 int cell_x
,cell_y
,cell_z
;
1773 int d
,tx
,ty
,tz
,dx
,dy
,dz
,cj
;
1774 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1775 int zsh_ty
,zsh_tx
,ysh_tx
;
1777 int dx0
,dx1
,dy0
,dy1
,dz0
,dz1
;
1778 int Nx
,Ny
,Nz
,shift
=-1,j
,nrj
,nns
,nn
=-1;
1779 real gridx
,gridy
,gridz
,grid_x
,grid_y
,grid_z
;
1780 real
*dcx2
,*dcy2
,*dcz2
;
1782 int cg0
,cg1
,icg
=-1,cgsnr
,i0
,igid
,nri
,naaj
,max_jcg
;
1783 int jcg0
,jcg1
,jjcg
,cgj0
,jgid
;
1784 int *grida
,*gridnra
,*gridind
;
1785 bool rvdw_lt_rcoul
,rcoul_lt_rvdw
;
1786 rvec xi
,*cgcm
,grid_offset
;
1787 real r2
,rs2
,rvdw2
,rcoul2
,rm2
,rl2
,XI
,YI
,ZI
,dcx
,dcy
,dcz
,tmp1
,tmp2
;
1789 bool bDomDec
,bTriclinicX
,bTriclinicY
;
1794 bDomDec
= DOMAINDECOMP(cr
);
1800 bTriclinicX
= ((YY
< grid
->npbcdim
&&
1801 (!bDomDec
|| dd
->nc
[YY
]==1) && box
[YY
][XX
] != 0) ||
1802 (ZZ
< grid
->npbcdim
&&
1803 (!bDomDec
|| dd
->nc
[ZZ
]==1) && box
[ZZ
][XX
] != 0));
1804 bTriclinicY
= (ZZ
< grid
->npbcdim
&&
1805 (!bDomDec
|| dd
->nc
[ZZ
]==1) && box
[ZZ
][YY
] != 0);
1808 rs2
= sqr(fr
->rlist
);
1809 if (bDoLongRange
&& fr
->bTwinRange
)
1811 /* The VdW and elec. LR cut-off's could be different,
1812 * so we can not simply set them to rlistlong.
1814 if (EVDW_ZERO_AT_CUTOFF(fr
->vdwtype
) && fr
->rvdw
> fr
->rlist
)
1816 rvdw2
= sqr(fr
->rlistlong
);
1820 rvdw2
= sqr(fr
->rvdw
);
1822 if (EEL_ZERO_AT_CUTOFF(fr
->eeltype
) && fr
->rcoulomb
> fr
->rlist
)
1824 rcoul2
= sqr(fr
->rlistlong
);
1828 rcoul2
= sqr(fr
->rcoulomb
);
1833 /* Workaround for a gcc -O3 or -ffast-math problem */
1837 rm2
= min(rvdw2
,rcoul2
);
1838 rl2
= max(rvdw2
,rcoul2
);
1839 rvdw_lt_rcoul
= (rvdw2
>= rcoul2
);
1840 rcoul_lt_rvdw
= (rcoul2
>= rvdw2
);
1842 if (bMakeQMMMnblist
)
1848 if (ns
->nl_sr
== NULL
)
1850 /* Short range buffers */
1851 snew(ns
->nl_sr
,ngid
);
1854 snew(ns
->nlr_ljc
,ngid
);
1855 snew(ns
->nlr_one
,ngid
);
1859 /* Long range VdW and Coul buffers */
1860 snew(ns
->nl_lr_ljc
,ngid
);
1864 /* Long range VdW or Coul only buffers */
1865 snew(ns
->nl_lr_one
,ngid
);
1867 for(j
=0; (j
<ngid
); j
++) {
1868 snew(ns
->nl_sr
[j
],MAX_CG
);
1871 snew(ns
->nl_lr_ljc
[j
],MAX_CG
);
1875 snew(ns
->nl_lr_one
[j
],MAX_CG
);
1880 "ns5_core: rs2 = %g, rvdw2 = %g, rcoul2 = %g (nm^2)\n",
1885 nl_lr_ljc
= ns
->nl_lr_ljc
;
1886 nl_lr_one
= ns
->nl_lr_one
;
1887 nlr_ljc
= ns
->nlr_ljc
;
1888 nlr_one
= ns
->nlr_one
;
1896 gridind
= grid
->index
;
1897 gridnra
= grid
->nra
;
1900 gridx
= grid
->cell_size
[XX
];
1901 gridy
= grid
->cell_size
[YY
];
1902 gridz
= grid
->cell_size
[ZZ
];
1906 copy_rvec(grid
->cell_offset
,grid_offset
);
1907 copy_ivec(grid
->ncpddc
,ncpddc
);
1912 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1913 zsh_ty
= floor(-box
[ZZ
][YY
]/box
[YY
][YY
]+0.5);
1914 zsh_tx
= floor(-box
[ZZ
][XX
]/box
[XX
][XX
]+0.5);
1915 ysh_tx
= floor(-box
[YY
][XX
]/box
[XX
][XX
]+0.5);
1916 if (zsh_tx
!=0 && ysh_tx
!=0)
1918 /* This could happen due to rounding, when both ratios are 0.5 */
1927 /* We only want a list for the test particle */
1936 /* Set the shift range */
1937 for(d
=0; d
<DIM
; d
++)
1941 /* Check if we need periodicity shifts.
1942 * Without PBC or with domain decomposition we don't need them.
1944 if (d
>= ePBC2npbcdim(fr
->ePBC
) || (bDomDec
&& dd
->nc
[d
] > 1))
1951 box
[XX
][XX
] - fabs(box
[YY
][XX
]) - fabs(box
[ZZ
][XX
]) < sqrt(rl2
))
1962 /* Loop over charge groups */
1963 for(icg
=cg0
; (icg
< cg1
); icg
++)
1965 igid
= GET_CGINFO_GID(cginfo
[icg
]);
1966 /* Skip this charge group if all energy groups are excluded! */
1967 if (bExcludeAlleg
[igid
])
1972 i0
= cgs
->index
[icg
];
1974 if (bMakeQMMMnblist
)
1976 /* Skip this charge group if it is not a QM atom while making a
1977 * QM/MM neighbourlist
1979 if (md
->bQM
[i0
]==FALSE
)
1981 continue; /* MM particle, go to next particle */
1984 /* Compute the number of charge groups that fall within the control
1987 naaj
= calc_naaj(icg
,cgsnr
);
1994 /* make a normal neighbourlist */
1998 /* Get the j charge-group and dd cell shift ranges */
1999 dd_get_ns_ranges(cr
->dd
,icg
,&jcg0
,&jcg1
,sh0
,sh1
);
2004 /* Compute the number of charge groups that fall within the control
2007 naaj
= calc_naaj(icg
,cgsnr
);
2013 /* The i-particle is awlways the test particle,
2014 * so we want all j-particles
2016 max_jcg
= cgsnr
- 1;
2020 max_jcg
= jcg1
- cgsnr
;
2025 i_egp_flags
= fr
->egp_flags
+ igid
*ngid
;
2027 /* Set the exclusions for the atoms in charge group icg using a bitmask */
2028 setexcl(i0
,cgs
->index
[icg
+1],&top
->excls
,TRUE
,bexcl
);
2030 ci2xyz(grid
,icg
,&cell_x
,&cell_y
,&cell_z
);
2032 /* Changed iicg to icg, DvdS 990115
2033 * (but see consistency check above, DvdS 990330)
2036 fprintf(log
,"icg=%5d, naaj=%5d, cell %d %d %d\n",
2037 icg
,naaj
,cell_x
,cell_y
,cell_z
);
2039 /* Loop over shift vectors in three dimensions */
2040 for (tz
=-shp
[ZZ
]; tz
<=shp
[ZZ
]; tz
++)
2042 ZI
= cgcm
[icg
][ZZ
]+tz
*box
[ZZ
][ZZ
];
2043 /* Calculate range of cells in Z direction that have the shift tz */
2044 zgi
= cell_z
+ tz
*Nz
;
2047 get_dx(Nz
,gridz
,rl2
,zgi
,ZI
,&dz0
,&dz1
,dcz2
);
2049 get_dx_dd(Nz
,gridz
,rl2
,zgi
,ZI
-grid_offset
[ZZ
],
2050 ncpddc
[ZZ
],sh0
[ZZ
],sh1
[ZZ
],&dz0
,&dz1
,dcz2
);
2056 for (ty
=-shp
[YY
]; ty
<=shp
[YY
]; ty
++)
2058 YI
= cgcm
[icg
][YY
]+ty
*box
[YY
][YY
]+tz
*box
[ZZ
][YY
];
2059 /* Calculate range of cells in Y direction that have the shift ty */
2062 ygi
= (int)(Ny
+ (YI
- grid_offset
[YY
])*grid_y
) - Ny
;
2066 ygi
= cell_y
+ ty
*Ny
;
2069 get_dx(Ny
,gridy
,rl2
,ygi
,YI
,&dy0
,&dy1
,dcy2
);
2071 get_dx_dd(Ny
,gridy
,rl2
,ygi
,YI
-grid_offset
[YY
],
2072 ncpddc
[YY
],sh0
[YY
],sh1
[YY
],&dy0
,&dy1
,dcy2
);
2078 for (tx
=-shp
[XX
]; tx
<=shp
[XX
]; tx
++)
2080 XI
= cgcm
[icg
][XX
]+tx
*box
[XX
][XX
]+ty
*box
[YY
][XX
]+tz
*box
[ZZ
][XX
];
2081 /* Calculate range of cells in X direction that have the shift tx */
2084 xgi
= (int)(Nx
+ (XI
- grid_offset
[XX
])*grid_x
) - Nx
;
2088 xgi
= cell_x
+ tx
*Nx
;
2091 get_dx(Nx
,gridx
,rl2
,xgi
*Nx
,XI
,&dx0
,&dx1
,dcx2
);
2093 get_dx_dd(Nx
,gridx
,rl2
,xgi
,XI
-grid_offset
[XX
],
2094 ncpddc
[XX
],sh0
[XX
],sh1
[XX
],&dx0
,&dx1
,dcx2
);
2100 /* Get shift vector */
2101 shift
=XYZ2IS(tx
,ty
,tz
);
2103 range_check(shift
,0,SHIFTS
);
2105 for(nn
=0; (nn
<ngid
); nn
++)
2112 fprintf(log
,"shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
2113 shift
,dx0
,dx1
,dy0
,dy1
,dz0
,dz1
);
2114 fprintf(log
,"cgcm: %8.3f %8.3f %8.3f\n",cgcm
[icg
][XX
],
2115 cgcm
[icg
][YY
],cgcm
[icg
][ZZ
]);
2116 fprintf(log
,"xi: %8.3f %8.3f %8.3f\n",XI
,YI
,ZI
);
2118 for (dx
=dx0
; (dx
<=dx1
); dx
++)
2120 tmp1
= rl2
- dcx2
[dx
];
2121 for (dy
=dy0
; (dy
<=dy1
); dy
++)
2123 tmp2
= tmp1
- dcy2
[dy
];
2126 for (dz
=dz0
; (dz
<=dz1
); dz
++) {
2127 if (tmp2
> dcz2
[dz
]) {
2128 /* Find grid-cell cj in which possible neighbours are */
2129 cj
= xyz2ci(Ny
,Nz
,dx
,dy
,dz
);
2131 /* Check out how many cgs (nrj) there in this cell */
2134 /* Find the offset in the cg list */
2137 /* Check if all j's are out of range so we
2138 * can skip the whole cell.
2139 * Should save some time, especially with DD.
2142 (grida
[cgj0
] >= max_jcg
&&
2143 (grida
[cgj0
] >= jcg1
|| grida
[cgj0
+nrj
-1] < jcg0
)))
2149 for (j
=0; (j
<nrj
); j
++)
2151 jjcg
= grida
[cgj0
+j
];
2153 /* check whether this guy is in range! */
2154 if ((jjcg
>= jcg0
&& jjcg
< jcg1
) ||
2157 r2
=calc_dx2(XI
,YI
,ZI
,cgcm
[jjcg
]);
2159 /* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
2160 jgid
= GET_CGINFO_GID(cginfo
[jjcg
]);
2161 /* check energy group exclusions */
2162 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
2166 if (nsr
[jgid
] >= MAX_CG
)
2168 put_in_list(bHaveVdW
,ngid
,md
,icg
,jgid
,
2169 nsr
[jgid
],nl_sr
[jgid
],
2170 cgs
->index
,/* cgsatoms, */ bexcl
,
2171 shift
,fr
,FALSE
,TRUE
,TRUE
,
2175 nl_sr
[jgid
][nsr
[jgid
]++]=jjcg
;
2179 if (nlr_ljc
[jgid
] >= MAX_CG
)
2181 do_longrange(cr
,top
,fr
,ngid
,md
,icg
,jgid
,
2183 nl_lr_ljc
[jgid
],bexcl
,shift
,x
,
2192 nl_lr_ljc
[jgid
][nlr_ljc
[jgid
]++]=jjcg
;
2196 if (nlr_one
[jgid
] >= MAX_CG
) {
2197 do_longrange(cr
,top
,fr
,ngid
,md
,icg
,jgid
,
2199 nl_lr_one
[jgid
],bexcl
,shift
,x
,
2203 rvdw_lt_rcoul
,rcoul_lt_rvdw
,FALSE
,
2208 nl_lr_one
[jgid
][nlr_one
[jgid
]++]=jjcg
;
2220 /* CHECK whether there is anything left in the buffers */
2221 for(nn
=0; (nn
<ngid
); nn
++)
2225 put_in_list(bHaveVdW
,ngid
,md
,icg
,nn
,nsr
[nn
],nl_sr
[nn
],
2226 cgs
->index
, /* cgsatoms, */ bexcl
,
2227 shift
,fr
,FALSE
,TRUE
,TRUE
,bMakeQMMMnblist
);
2230 if (nlr_ljc
[nn
] > 0)
2232 do_longrange(cr
,top
,fr
,ngid
,md
,icg
,nn
,nlr_ljc
[nn
],
2233 nl_lr_ljc
[nn
],bexcl
,shift
,x
,box_size
,nrnb
,
2234 lambda
,dvdlambda
,grppener
,TRUE
,TRUE
,FALSE
,
2235 bHaveVdW
,bDoForces
,f
);
2238 if (nlr_one
[nn
] > 0)
2240 do_longrange(cr
,top
,fr
,ngid
,md
,icg
,nn
,nlr_one
[nn
],
2241 nl_lr_one
[nn
],bexcl
,shift
,x
,box_size
,nrnb
,
2242 lambda
,dvdlambda
,grppener
,
2243 rvdw_lt_rcoul
,rcoul_lt_rvdw
,FALSE
,
2244 bHaveVdW
,bDoForces
,f
);
2250 /* setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl); */
2251 setexcl(cgs
->index
[icg
],cgs
->index
[icg
+1],&top
->excls
,FALSE
,bexcl
);
2253 /* Perform any left over force calculations */
2254 for (nn
=0; (nn
<ngid
); nn
++)
2258 do_longrange(cr
,top
,fr
,0,md
,icg
,nn
,nlr_ljc
[nn
],
2259 nl_lr_ljc
[nn
],bexcl
,shift
,x
,box_size
,nrnb
,
2260 lambda
,dvdlambda
,grppener
,
2261 TRUE
,TRUE
,TRUE
,bHaveVdW
,bDoForces
,f
);
2264 do_longrange(cr
,top
,fr
,0,md
,icg
,nn
,nlr_one
[nn
],
2265 nl_lr_one
[nn
],bexcl
,shift
,x
,box_size
,nrnb
,
2266 lambda
,dvdlambda
,grppener
,
2267 rvdw_lt_rcoul
,rcoul_lt_rvdw
,
2268 TRUE
,bHaveVdW
,bDoForces
,f
);
2273 /* Close off short range neighbourlists */
2274 close_neighbor_list(fr
,FALSE
,-1,-1,bMakeQMMMnblist
);
2279 void ns_realloc_natoms(gmx_ns_t
*ns
,int natoms
)
2283 if (natoms
> ns
->nra_alloc
)
2285 ns
->nra_alloc
= over_alloc_dd(natoms
);
2286 srenew(ns
->bexcl
,ns
->nra_alloc
);
2287 for(i
=0; i
<ns
->nra_alloc
; i
++)
2294 void init_ns(FILE *fplog
,const t_commrec
*cr
,
2295 gmx_ns_t
*ns
,t_forcerec
*fr
,
2296 const gmx_mtop_t
*mtop
,
2299 int mt
,icg
,nr_in_cg
,maxcg
,i
,j
,jcg
,ngid
,ncg
;
2303 /* Compute largest charge groups size (# atoms) */
2305 for(mt
=0; mt
<mtop
->nmoltype
; mt
++) {
2306 cgs
= &mtop
->moltype
[mt
].cgs
;
2307 for (icg
=0; (icg
< cgs
->nr
); icg
++)
2309 nr_in_cg
=max(nr_in_cg
,(int)(cgs
->index
[icg
+1]-cgs
->index
[icg
]));
2313 /* Verify whether largest charge group is <= max cg.
2314 * This is determined by the type of the local exclusion type
2315 * Exclusions are stored in bits. (If the type is not large
2316 * enough, enlarge it, unsigned char -> unsigned short -> unsigned long)
2318 maxcg
= sizeof(t_excl
)*8;
2319 if (nr_in_cg
> maxcg
)
2321 gmx_fatal(FARGS
,"Max #atoms in a charge group: %d > %d\n",
2325 ngid
= mtop
->groups
.grps
[egcENER
].nr
;
2326 snew(ns
->bExcludeAlleg
,ngid
);
2327 for(i
=0; i
<ngid
; i
++) {
2328 ns
->bExcludeAlleg
[i
] = TRUE
;
2329 for(j
=0; j
<ngid
; j
++)
2331 if (!(fr
->egp_flags
[i
*ngid
+j
] & EGP_EXCL
))
2333 ns
->bExcludeAlleg
[i
] = FALSE
;
2340 ns
->grid
= init_grid(fplog
,fr
);
2341 /* These lists are allocated in ns5_core */
2347 snew(ns
->ns_buf
,ngid
);
2348 for(i
=0; (i
<ngid
); i
++)
2350 snew(ns
->ns_buf
[i
],SHIFTS
);
2352 ncg
= ncg_mtop(mtop
);
2353 snew(ns
->simple_aaj
,2*ncg
);
2354 for(jcg
=0; (jcg
<ncg
); jcg
++)
2356 ns
->simple_aaj
[jcg
] = jcg
;
2357 ns
->simple_aaj
[jcg
+ncg
] = jcg
;
2361 /* Create array that determines whether or not atoms have VdW */
2362 snew(ns
->bHaveVdW
,fr
->ntype
);
2363 for(i
=0; (i
<fr
->ntype
); i
++)
2365 for(j
=0; (j
<fr
->ntype
); j
++)
2367 ns
->bHaveVdW
[i
] = (ns
->bHaveVdW
[i
] ||
2369 ((BHAMA(fr
->nbfp
,fr
->ntype
,i
,j
) != 0) ||
2370 (BHAMB(fr
->nbfp
,fr
->ntype
,i
,j
) != 0) ||
2371 (BHAMC(fr
->nbfp
,fr
->ntype
,i
,j
) != 0)) :
2372 ((C6(fr
->nbfp
,fr
->ntype
,i
,j
) != 0) ||
2373 (C12(fr
->nbfp
,fr
->ntype
,i
,j
) != 0))));
2377 pr_bvec(debug
,0,"bHaveVdW",ns
->bHaveVdW
,fr
->ntype
,TRUE
);
2379 if (!DOMAINDECOMP(cr
))
2381 /* This could be reduced with particle decomposition */
2382 ns_realloc_natoms(ns
,mtop
->natoms
);
2385 ns
->nblist_initialized
=FALSE
;
2387 /* nbr list debug dump */
2389 char *ptr
=getenv("GMX_DUMP_NL");
2392 ns
->dump_nl
=strtol(ptr
,NULL
,0);
2395 fprintf(fplog
, "GMX_DUMP_NL = %d", ns
->dump_nl
);
2406 int search_neighbours(FILE *log
,t_forcerec
*fr
,
2407 rvec x
[],matrix box
,
2408 gmx_localtop_t
*top
,
2409 gmx_groups_t
*groups
,
2411 t_nrnb
*nrnb
,t_mdatoms
*md
,
2412 real lambda
,real
*dvdlambda
,
2413 gmx_grppairener_t
*grppener
,
2416 bool bDoForces
,rvec
*f
)
2418 t_block
*cgs
=&(top
->cgs
);
2419 rvec box_size
,grid_x0
,grid_x1
;
2421 real min_size
,grid_dens
;
2426 int cg_start
,cg_end
,start
,end
;
2429 gmx_domdec_zones_t
*dd_zones
;
2433 /* Set some local variables */
2435 ngid
= groups
->grps
[egcENER
].nr
;
2437 for(m
=0; (m
<DIM
); m
++)
2439 box_size
[m
] = box
[m
][m
];
2442 if (fr
->ePBC
!= epbcNONE
)
2444 if (sqr(fr
->rlistlong
) >= max_cutoff2(fr
->ePBC
,box
))
2446 gmx_fatal(FARGS
,"One of the box vectors has become shorter than twice the cut-off length or box_yy-|box_zy| or box_zz has become smaller than the cut-off.");
2450 min_size
= min(box_size
[XX
],min(box_size
[YY
],box_size
[ZZ
]));
2451 if (2*fr
->rlistlong
>= min_size
)
2452 gmx_fatal(FARGS
,"One of the box diagonal elements has become smaller than twice the cut-off length.");
2456 if (DOMAINDECOMP(cr
))
2458 ns_realloc_natoms(ns
,cgs
->index
[cgs
->nr
]);
2462 /* Reset the neighbourlists */
2463 reset_neighbor_list(fr
,FALSE
,-1,-1);
2465 if (bGrid
&& bFillGrid
)
2469 if (DOMAINDECOMP(cr
))
2471 dd_zones
= domdec_zones(cr
->dd
);
2477 get_nsgrid_boundaries(grid
,NULL
,box
,NULL
,NULL
,NULL
,
2478 cgs
->nr
,fr
->cg_cm
,grid_x0
,grid_x1
,&grid_dens
);
2480 grid_first(log
,grid
,NULL
,NULL
,fr
->ePBC
,box
,grid_x0
,grid_x1
,
2481 fr
->rlistlong
,grid_dens
);
2485 /* Don't know why this all is... (DvdS 3/99) */
2491 end
= (cgs
->nr
+1)/2;
2494 if (DOMAINDECOMP(cr
))
2497 fill_grid(log
,dd_zones
,grid
,end
,-1,end
,fr
->cg_cm
);
2499 grid
->icg1
= dd_zones
->izone
[dd_zones
->nizone
-1].cg1
;
2503 fill_grid(log
,NULL
,grid
,cgs
->nr
,fr
->cg0
,fr
->hcg
,fr
->cg_cm
);
2504 grid
->icg0
= fr
->cg0
;
2505 grid
->icg1
= fr
->hcg
;
2513 calc_elemnr(log
,grid
,start
,end
,cgs
->nr
);
2515 grid_last(log
,grid
,start
,end
,cgs
->nr
);
2519 check_grid(debug
,grid
);
2520 print_grid(debug
,grid
);
2525 /* Set the grid cell index for the test particle only.
2526 * The cell to cg index is not corrected, but that does not matter.
2528 fill_grid(log
,NULL
,ns
->grid
,fr
->hcg
,fr
->hcg
-1,fr
->hcg
,fr
->cg_cm
);
2536 nsearch
= nsgrid_core(log
,cr
,fr
,box
,box_size
,ngid
,top
,
2537 grid
,x
,ns
->bexcl
,ns
->bExcludeAlleg
,
2538 nrnb
,md
,lambda
,dvdlambda
,grppener
,
2540 bDoLongRange
,bDoForces
,f
,
2543 /* neighbour searching withouth QMMM! QM atoms have zero charge in
2544 * the classical calculation. The charge-charge interaction
2545 * between QM and MM atoms is handled in the QMMM core calculation
2546 * (see QMMM.c). The VDW however, we'd like to compute classically
2547 * and the QM MM atom pairs have just been put in the
2548 * corresponding neighbourlists. in case of QMMM we still need to
2549 * fill a special QMMM neighbourlist that contains all neighbours
2550 * of the QM atoms. If bQMMM is true, this list will now be made:
2552 if (fr
->bQMMM
&& fr
->qr
->QMMMscheme
!=eQMMMschemeoniom
)
2554 nsearch
+= nsgrid_core(log
,cr
,fr
,box
,box_size
,ngid
,top
,
2555 grid
,x
,ns
->bexcl
,ns
->bExcludeAlleg
,
2556 nrnb
,md
,lambda
,dvdlambda
,grppener
,
2558 bDoLongRange
,bDoForces
,f
,
2564 nsearch
= ns_simple_core(fr
,top
,md
,box
,box_size
,
2565 ns
->bexcl
,ns
->simple_aaj
,
2566 ngid
,ns
->ns_buf
,ns
->bHaveVdW
);
2574 inc_nrnb(nrnb
,eNR_NS
,nsearch
);
2575 /* inc_nrnb(nrnb,eNR_LR,fr->nlr); */
2580 int natoms_beyond_ns_buffer(t_inputrec
*ir
,t_forcerec
*fr
,t_block
*cgs
,
2581 matrix scale_tot
,rvec
*x
)
2583 int cg0
,cg1
,cg
,a0
,a1
,a
,i
,j
;
2584 real rint
,hbuf2
,scale
;
2591 rint
= max(ir
->rcoulomb
,ir
->rvdw
);
2592 if (ir
->rlist
< rint
)
2594 gmx_fatal(FARGS
,"The neighbor search buffer has negative size: %f nm",
2602 if (!EI_DYNAMICS(ir
->eI
) || !DYNAMIC_BOX(*ir
))
2604 hbuf2
= sqr(0.5*(ir
->rlist
- rint
));
2605 for(cg
=cg0
; cg
<cg1
; cg
++)
2607 a0
= cgs
->index
[cg
];
2608 a1
= cgs
->index
[cg
+1];
2609 for(a
=a0
; a
<a1
; a
++)
2611 if (distance2(cg_cm
[cg
],x
[a
]) > hbuf2
)
2621 scale
= scale_tot
[0][0];
2622 for(i
=1; i
<DIM
; i
++)
2624 /* With anisotropic scaling, the original spherical ns volumes become
2625 * ellipsoids. To avoid costly transformations we use the minimum
2626 * eigenvalue of the scaling matrix for determining the buffer size.
2627 * Since the lower half is 0, the eigenvalues are the diagonal elements.
2629 scale
= min(scale
,scale_tot
[i
][i
]);
2630 if (scale_tot
[i
][i
] != scale_tot
[i
-1][i
-1])
2636 if (scale_tot
[i
][j
] != 0)
2642 hbuf2
= sqr(0.5*(scale
*ir
->rlist
- rint
));
2645 for(cg
=cg0
; cg
<cg1
; cg
++)
2647 svmul(scale
,cg_cm
[cg
],cgsc
);
2648 a0
= cgs
->index
[cg
];
2649 a1
= cgs
->index
[cg
+1];
2650 for(a
=a0
; a
<a1
; a
++)
2652 if (distance2(cgsc
,x
[a
]) > hbuf2
)
2661 /* Anistropic scaling */
2662 for(cg
=cg0
; cg
<cg1
; cg
++)
2664 /* Since scale_tot contains the transpose of the scaling matrix,
2665 * we need to multiply with the transpose.
2667 tmvmul_ur0(scale_tot
,cg_cm
[cg
],cgsc
);
2668 a0
= cgs
->index
[cg
];
2669 a1
= cgs
->index
[cg
+1];
2670 for(a
=a0
; a
<a1
; a
++)
2672 if (distance2(cgsc
,x
[a
]) > hbuf2
)