?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
minor fixes in ditribution files
[gromacs/qmmm-gamess-us.git] / src / mdlib / domdec_top.c
blobd3f974d067c91bee3a88ae4fc425b11981b37333
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This file is part of Gromacs Copyright (c) 1991-2008
5 * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
6 *
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
11 *
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
14 *
15 * And Hey:
16 * Gnomes, ROck Monsters And Chili Sauce
17 */
18
19 #ifdef HAVE_CONFIG_H
20 #include <config.h>
21 #endif
22
23 #include <string.h>
24 #include "typedefs.h"
25 #include "smalloc.h"
26 #include "domdec.h"
27 #include "domdec_network.h"
28 #include "names.h"
29 #include "network.h"
30 #include "vec.h"
31 #include "pbc.h"
32 #include "chargegroup.h"
33 #include "gmx_random.h"
34 #include "topsort.h"
35 #include "mtop_util.h"
36 #include "mshift.h"
37 #include "vsite.h"
38 #include "gmx_ga2la.h"
39
40 typedef struct {
41 int *index; /* Index for each atom into il */
42 int *il; /* ftype|type|a0|...|an|ftype|... */
43 } gmx_reverse_ilist_t;
44
45 typedef struct {
46 int a_start;
47 int a_end;
48 int natoms_mol;
49 int type;
50 } gmx_molblock_ind_t;
51
52 typedef struct gmx_reverse_top {
53 bool bExclRequired; /* Do we require all exclusions to be assigned? */
54 bool bConstr; /* Are there constraints in this revserse top? */
55 bool bBCheck; /* All bonded interactions have to be assigned? */
56 bool bMultiCGmols; /* Are the multi charge-group molecules? */
57 gmx_reverse_ilist_t *ril_mt; /* Reverse ilist for all moltypes */
58 int ril_mt_tot_size;
59 int ilsort; /* The sorting state of bondeds for free energy */
60 gmx_molblock_ind_t *mbi;
61
62 /* Pointers only used for an error message */
63 gmx_mtop_t *err_top_global;
64 gmx_localtop_t *err_top_local;
65 } gmx_reverse_top_t;
66
67 static int nral_rt(int ftype)
68 {
69 /* Returns the number of atom entries for il in gmx_reverse_top_t */
70 int nral;
71
72 nral = NRAL(ftype);
73 if (interaction_function[ftype].flags & IF_VSITE)
74 {
75 /* With vsites the reverse topology contains
76 * two extra entries for PBC.
77 */
78 nral += 2;
79 }
80
81 return nral;
82 }
83
84 static bool dd_check_ftype(int ftype,bool bBCheck,bool bConstr)
85 {
86 return (((interaction_function[ftype].flags & IF_BOND) &&
87 !(interaction_function[ftype].flags & IF_VSITE) &&
88 (bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO))) ||
89 ftype == F_SETTLE ||
90 (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)));
91 }
92
93 static void print_error_header(FILE *fplog,char *moltypename,int nprint)
94 {
95 fprintf(fplog, "\nMolecule type '%s'\n",moltypename);
96 fprintf(stderr,"\nMolecule type '%s'\n",moltypename);
97 fprintf(fplog,
98 "the first %d missing interactions, except for exclusions:\n",
99 nprint);
100 fprintf(stderr,
101 "the first %d missing interactions, except for exclusions:\n",
102 nprint);
103 }
104
105 static void print_missing_interactions_mb(FILE *fplog,t_commrec *cr,
106 gmx_reverse_top_t *rt,
107 char *moltypename,
108 gmx_reverse_ilist_t *ril,
109 int a_start,int a_end,
110 int nat_mol,int nmol,
111 t_idef *idef)
112 {
113 int nril_mol,*assigned,*gatindex;
114 int ftype,ftype_j,nral,i,j_mol,j,k,a0,a0_mol,mol,a,a_gl;
115 int nprint;
116 t_ilist *il;
117 t_iatom *ia;
118 bool bFound;
119
120 nril_mol = ril->index[nat_mol];
121 snew(assigned,nmol*nril_mol);
122
123 gatindex = cr->dd->gatindex;
124 for(ftype=0; ftype<F_NRE; ftype++)
125 {
126 if (dd_check_ftype(ftype,rt->bBCheck,rt->bConstr))
127 {
128 nral = NRAL(ftype);
129 il = &idef->il[ftype];
130 ia = il->iatoms;
131 for(i=0; i<il->nr; i+=1+nral)
132 {
133 a0 = gatindex[ia[1]];
134 /* Check if this interaction is in
135 * the currently checked molblock.
136 */
137 if (a0 >= a_start && a0 < a_end)
138 {
139 mol = (a0 - a_start)/nat_mol;
140 a0_mol = (a0 - a_start) - mol*nat_mol;
141 j_mol = ril->index[a0_mol];
142 bFound = FALSE;
143 while (j_mol < ril->index[a0_mol+1] && !bFound)
144 {
145 j = mol*nril_mol + j_mol;
146 ftype_j = ril->il[j_mol];
147 /* Here we need to check if this interaction has
148 * not already been assigned, since we could have
149 * multiply defined interactions.
150 */
151 if (ftype == ftype_j && ia[0] == ril->il[j_mol+1] &&
152 assigned[j] == 0)
153 {
154 /* Check the atoms */
155 bFound = TRUE;
156 for(a=0; a<nral; a++)
157 {
158 if (gatindex[ia[1+a]] !=
159 a_start + mol*nat_mol + ril->il[j_mol+2+a])
160 {
161 bFound = FALSE;
162 }
163 }
164 if (bFound)
165 {
166 assigned[j] = 1;
167 }
168 }
169 j_mol += 2 + nral_rt(ftype_j);
170 }
171 if (!bFound)
172 {
173 gmx_incons("Some interactions seem to be assigned multiple times");
174 }
175 }
176 ia += 1 + nral;
177 }
178 }
179 }
180
181 gmx_sumi(nmol*nril_mol,assigned,cr);
182
183 nprint = 10;
184 i = 0;
185 for(mol=0; mol<nmol; mol++)
186 {
187 j_mol = 0;
188 while (j_mol < nril_mol)
189 {
190 ftype = ril->il[j_mol];
191 nral = NRAL(ftype);
192 j = mol*nril_mol + j_mol;
193 if (assigned[j] == 0 &&
194 !(interaction_function[ftype].flags & IF_VSITE))
195 {
196 if (DDMASTER(cr->dd))
197 {
198 if (i == 0)
199 {
200 print_error_header(fplog,moltypename,nprint);
201 }
202 fprintf(fplog, "%20s atoms",
203 interaction_function[ftype].longname);
204 fprintf(stderr,"%20s atoms",
205 interaction_function[ftype].longname);
206 for(a=0; a<nral; a++) {
207 fprintf(fplog, "%5d",ril->il[j_mol+2+a]+1);
208 fprintf(stderr,"%5d",ril->il[j_mol+2+a]+1);
209 }
210 while (a < 4)
211 {
212 fprintf(fplog, " ");
213 fprintf(stderr," ");
214 a++;
215 }
216 fprintf(fplog, " global");
217 fprintf(stderr," global");
218 for(a=0; a<nral; a++)
219 {
220 fprintf(fplog, "%6d",
221 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
222 fprintf(stderr,"%6d",
223 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
224 }
225 fprintf(fplog, "\n");
226 fprintf(stderr,"\n");
227 }
228 i++;
229 if (i >= nprint)
230 {
231 break;
232 }
233 }
234 j_mol += 2 + nral_rt(ftype);
235 }
236 }
237
238 sfree(assigned);
239 }
240
241 static void print_missing_interactions_atoms(FILE *fplog,t_commrec *cr,
242 gmx_mtop_t *mtop,t_idef *idef)
243 {
244 int mb,a_start,a_end;
245 gmx_molblock_t *molb;
246 gmx_reverse_top_t *rt;
247
248 rt = cr->dd->reverse_top;
249
250 /* Print the atoms in the missing interactions per molblock */
251 a_end = 0;
252 for(mb=0; mb<mtop->nmolblock; mb++)
253 {
254 molb = &mtop->molblock[mb];
255 a_start = a_end;
256 a_end = a_start + molb->nmol*molb->natoms_mol;
257
258 print_missing_interactions_mb(fplog,cr,rt,
259 *(mtop->moltype[molb->type].name),
260 &rt->ril_mt[molb->type],
261 a_start,a_end,molb->natoms_mol,
262 molb->nmol,
263 idef);
264 }
265 }
266
267 void dd_print_missing_interactions(FILE *fplog,t_commrec *cr,int local_count, gmx_mtop_t *top_global, t_state *state_local)
268 {
269 int ndiff_tot,cl[F_NRE],n,ndiff,rest_global,rest_local;
270 int ftype,nral;
271 char buf[STRLEN];
272 gmx_domdec_t *dd;
273 gmx_mtop_t *err_top_global;
274 gmx_localtop_t *err_top_local;
275
276 dd = cr->dd;
277
278 err_top_global = dd->reverse_top->err_top_global;
279 err_top_local = dd->reverse_top->err_top_local;
280
281 if (fplog)
282 {
283 fprintf(fplog,"\nNot all bonded interactions have been properly assigned to the domain decomposition cells\n");
284 fflush(fplog);
285 }
286
287 ndiff_tot = local_count - dd->nbonded_global;
288
289 for(ftype=0; ftype<F_NRE; ftype++)
290 {
291 nral = NRAL(ftype);
292 cl[ftype] = err_top_local->idef.il[ftype].nr/(1+nral);
293 }
294
295 gmx_sumi(F_NRE,cl,cr);
296
297 if (DDMASTER(dd))
298 {
299 fprintf(fplog,"\nA list of missing interactions:\n");
300 fprintf(stderr,"\nA list of missing interactions:\n");
301 rest_global = dd->nbonded_global;
302 rest_local = local_count;
303 for(ftype=0; ftype<F_NRE; ftype++)
304 {
305 /* In the reverse and local top all constraints are merged
306 * into F_CONSTR. So in the if statement we skip F_CONSTRNC
307 * and add these constraints when doing F_CONSTR.
308 */
309 if (((interaction_function[ftype].flags & IF_BOND) &&
310 (dd->reverse_top->bBCheck
311 || !(interaction_function[ftype].flags & IF_LIMZERO)))
312 || ftype == F_SETTLE
313 || (dd->reverse_top->bConstr && ftype == F_CONSTR))
314 {
315 nral = NRAL(ftype);
316 n = gmx_mtop_ftype_count(err_top_global,ftype);
317 if (ftype == F_CONSTR)
318 {
319 n += gmx_mtop_ftype_count(err_top_global,F_CONSTRNC);
320 }
321 ndiff = cl[ftype] - n;
322 if (ndiff != 0)
323 {
324 sprintf(buf,"%20s of %6d missing %6d",
325 interaction_function[ftype].longname,n,-ndiff);
326 fprintf(fplog,"%s\n",buf);
327 fprintf(stderr,"%s\n",buf);
328 }
329 rest_global -= n;
330 rest_local -= cl[ftype];
331 }
332 }
333
334 ndiff = rest_local - rest_global;
335 if (ndiff != 0)
336 {
337 sprintf(buf,"%20s of %6d missing %6d","exclusions",
338 rest_global,-ndiff);
339 fprintf(fplog,"%s\n",buf);
340 fprintf(stderr,"%s\n",buf);
341 }
342 }
343
344 print_missing_interactions_atoms(fplog,cr,err_top_global,
345 &err_top_local->idef);
346 write_dd_pdb("dd_dump_err",0,"dump",top_global,cr,
347 -1,state_local->x,state_local->box);
348 if (DDMASTER(dd))
349 {
350 if (ndiff_tot > 0)
351 {
352 gmx_incons("One or more interactions were multiple assigned in the domain decompostion");
353 }
354 else
355 {
356 gmx_fatal(FARGS,"%d of the %d bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (%g nm) or the two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck",-ndiff_tot,cr->dd->nbonded_global,dd_cutoff_mbody(cr->dd),dd_cutoff_twobody(cr->dd));
357 }
358 }
359 }
360
361 static void global_atomnr_to_moltype_ind(gmx_molblock_ind_t *mbi,int i_gl,
362 int *mb,int *mt,int *mol,int *i_mol)
363 {
364 int molb;
365
366 *mb = 0;
367 while (i_gl >= mbi->a_end) {
368 (*mb)++;
369 mbi++;
370 }
371
372 *mt = mbi->type;
373 *mol = (i_gl - mbi->a_start) / mbi->natoms_mol;
374 *i_mol = (i_gl - mbi->a_start) - (*mol)*mbi->natoms_mol;
375 }
376
377 static int count_excls(t_block *cgs,t_blocka *excls,int *n_intercg_excl)
378 {
379 int n,n_inter,cg,at0,at1,at,excl,atj;
380
381 n = 0;
382 *n_intercg_excl = 0;
383 for(cg=0; cg<cgs->nr; cg++)
384 {
385 at0 = cgs->index[cg];
386 at1 = cgs->index[cg+1];
387 for(at=at0; at<at1; at++)
388 {
389 for(excl=excls->index[at]; excl<excls->index[at+1]; excl++)
390 {
391 atj = excls->a[excl];
392 if (atj > at)
393 {
394 n++;
395 if (atj < at0 || atj >= at1)
396 {
397 (*n_intercg_excl)++;
398 }
399 }
400 }
401 }
402 }
403
404 return n;
405 }
406
407 static int low_make_reverse_ilist(t_ilist *il_mt,t_atom *atom,
408 int **vsite_pbc,
409 int *count,
410 bool bConstr,bool bBCheck,
411 int *r_index,int *r_il,
412 bool bLinkToAllAtoms,
413 bool bAssign)
414 {
415 int ftype,nral,i,j,nlink,link;
416 t_ilist *il;
417 t_iatom *ia;
418 atom_id a;
419 int nint;
420 bool bVSite;
421
422 nint = 0;
423 for(ftype=0; ftype<F_NRE; ftype++)
424 {
425 if ((interaction_function[ftype].flags & (IF_BOND | IF_VSITE)) ||
426 ftype == F_SETTLE ||
427 (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC))) {
428 bVSite = (interaction_function[ftype].flags & IF_VSITE);
429 nral = NRAL(ftype);
430 il = &il_mt[ftype];
431 ia = il->iatoms;
432 for(i=0; i<il->nr; i+=1+nral)
433 {
434 ia = il->iatoms + i;
435 if (bLinkToAllAtoms)
436 {
437 if (bVSite)
438 {
439 /* We don't need the virtual sites for the cg-links */
440 nlink = 0;
441 }
442 else
443 {
444 nlink = nral;
445 }
446 }
447 else
448 {
449 /* Couple to the first atom in the interaction */
450 nlink = 1;
451 }
452 for(link=0; link<nlink; link++)
453 {
454 a = ia[1+link];
455 if (bAssign)
456 {
457 r_il[r_index[a]+count[a]] =
458 (ftype == F_CONSTRNC ? F_CONSTR : ftype);
459 r_il[r_index[a]+count[a]+1] = ia[0];
460 for(j=1; j<1+nral; j++)
461 {
462 /* Store the molecular atom number */
463 r_il[r_index[a]+count[a]+1+j] = ia[j];
464 }
465 }
466 if (interaction_function[ftype].flags & IF_VSITE)
467 {
468 if (bAssign)
469 {
470 /* Add an entry to iatoms for storing
471 * which of the constructing atoms are
472 * vsites again.
473 */
474 r_il[r_index[a]+count[a]+2+nral] = 0;
475 for(j=2; j<1+nral; j++)
476 {
477 if (atom[ia[j]].ptype == eptVSite)
478 {
479 r_il[r_index[a]+count[a]+2+nral] |= (2<<j);
480 }
481 }
482 /* Store vsite pbc atom in a second extra entry */
483 r_il[r_index[a]+count[a]+2+nral+1] =
484 (vsite_pbc ? vsite_pbc[ftype-F_VSITE2][i/(1+nral)] : -2);
485 }
486 }
487 else
488 {
489 /* We do not count vsites since they are always
490 * uniquely assigned and can be assigned
491 * to multiple nodes with recursive vsites.
492 */
493 if (bBCheck ||
494 !(interaction_function[ftype].flags & IF_LIMZERO))
495 {
496 nint++;
497 }
498 }
499 count[a] += 2 + nral_rt(ftype);
500 }
501 }
502 }
503 }
504
505 return nint;
506 }
507
508 static int make_reverse_ilist(gmx_moltype_t *molt,
509 int **vsite_pbc,
510 bool bConstr,bool bBCheck,
511 bool bLinkToAllAtoms,
512 gmx_reverse_ilist_t *ril_mt)
513 {
514 int nat_mt,*count,i,nint_mt;
515
516 /* Count the interactions */
517 nat_mt = molt->atoms.nr;
518 snew(count,nat_mt);
519 low_make_reverse_ilist(molt->ilist,molt->atoms.atom,vsite_pbc,
520 count,
521 bConstr,bBCheck,NULL,NULL,
522 bLinkToAllAtoms,FALSE);
523
524 snew(ril_mt->index,nat_mt+1);
525 ril_mt->index[0] = 0;
526 for(i=0; i<nat_mt; i++)
527 {
528 ril_mt->index[i+1] = ril_mt->index[i] + count[i];
529 count[i] = 0;
530 }
531 snew(ril_mt->il,ril_mt->index[nat_mt]);
532
533 /* Store the interactions */
534 nint_mt =
535 low_make_reverse_ilist(molt->ilist,molt->atoms.atom,vsite_pbc,
536 count,
537 bConstr,bBCheck,
538 ril_mt->index,ril_mt->il,
539 bLinkToAllAtoms,TRUE);
540
541 sfree(count);
542
543 return nint_mt;
544 }
545
546 static void destroy_reverse_ilist(gmx_reverse_ilist_t *ril)
547 {
548 sfree(ril->index);
549 sfree(ril->il);
550 }
551
552 static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop,bool bFE,
553 int ***vsite_pbc_molt,
554 bool bConstr,
555 bool bBCheck,int *nint)
556 {
557 int mt,i,mb;
558 gmx_reverse_top_t *rt;
559 int *nint_mt;
560 gmx_moltype_t *molt;
561
562 snew(rt,1);
563
564 /* Should we include constraints (for SHAKE) in rt? */
565 rt->bConstr = bConstr;
566 rt->bBCheck = bBCheck;
567
568 rt->bMultiCGmols = FALSE;
569 snew(nint_mt,mtop->nmoltype);
570 snew(rt->ril_mt,mtop->nmoltype);
571 rt->ril_mt_tot_size = 0;
572 for(mt=0; mt<mtop->nmoltype; mt++)
573 {
574 molt = &mtop->moltype[mt];
575 if (molt->cgs.nr > 1)
576 {
577 rt->bMultiCGmols = TRUE;
578 }
579
580 /* Make the atom to interaction list for this molecule type */
581 nint_mt[mt] =
582 make_reverse_ilist(molt,vsite_pbc_molt ? vsite_pbc_molt[mt] : NULL,
583 rt->bConstr,rt->bBCheck,FALSE,
584 &rt->ril_mt[mt]);
585
586 rt->ril_mt_tot_size += rt->ril_mt[mt].index[molt->atoms.nr];
587 }
588 if (debug)
589 {
590 fprintf(debug,"The total size of the atom to interaction index is %d integers\n",rt->ril_mt_tot_size);
591 }
592
593 *nint = 0;
594 for(mb=0; mb<mtop->nmolblock; mb++)
595 {
596 *nint += mtop->molblock[mb].nmol*nint_mt[mtop->molblock[mb].type];
597 }
598 sfree(nint_mt);
599
600 if (bFE && gmx_mtop_bondeds_free_energy(mtop))
601 {
602 rt->ilsort = ilsortFE_UNSORTED;
603 }
604 else {
605 rt->ilsort = ilsortNO_FE;
606 }
607
608 /* Make a molblock index for fast searching */
609 snew(rt->mbi,mtop->nmolblock);
610 i = 0;
611 for(mb=0; mb<mtop->nmolblock; mb++)
612 {
613 rt->mbi[mb].a_start = i;
614 i += mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
615 rt->mbi[mb].a_end = i;
616 rt->mbi[mb].natoms_mol = mtop->molblock[mb].natoms_mol;
617 rt->mbi[mb].type = mtop->molblock[mb].type;
618 }
619
620 return rt;
621 }
622
623 void dd_make_reverse_top(FILE *fplog,
624 gmx_domdec_t *dd,gmx_mtop_t *mtop,
625 gmx_vsite_t *vsite,gmx_constr_t constr,
626 t_inputrec *ir,bool bBCheck)
627 {
628 int mb,natoms,n_recursive_vsite,nexcl,nexcl_icg,a;
629 gmx_molblock_t *molb;
630 gmx_moltype_t *molt;
631
632 if (fplog)
633 {
634 fprintf(fplog,"\nLinking all bonded interactions to atoms\n");
635 }
636
637 dd->reverse_top = make_reverse_top(mtop,ir->efep!=efepNO,
638 vsite ? vsite->vsite_pbc_molt : NULL,
639 !dd->bInterCGcons,
640 bBCheck,&dd->nbonded_global);
641
642 if (dd->reverse_top->ril_mt_tot_size >= 200000 &&
643 mtop->mols.nr > 1 &&
644 mtop->nmolblock == 1 && mtop->molblock[0].nmol == 1)
645 {
646 /* mtop comes from a pre Gromacs 4 tpr file */
647 const char *note="NOTE: The tpr file used for this simulation is in an old format, for less memory usage and possibly more performance create a new tpr file with an up to date version of grompp";
648 if (fplog)
649 {
650 fprintf(fplog,"\n%s\n\n",note);
651 }
652 if (DDMASTER(dd))
653 {
654 fprintf(stderr,"\n%s\n\n",note);
655 }
656 }
657
658 dd->reverse_top->bExclRequired = IR_EXCL_FORCES(*ir);
659
660 nexcl = 0;
661 dd->n_intercg_excl = 0;
662 for(mb=0; mb<mtop->nmolblock; mb++)
663 {
664 molb = &mtop->molblock[mb];
665 molt = &mtop->moltype[molb->type];
666 nexcl += molb->nmol*count_excls(&molt->cgs,&molt->excls,&nexcl_icg);
667 dd->n_intercg_excl += molb->nmol*nexcl_icg;
668 }
669 if (dd->reverse_top->bExclRequired)
670 {
671 dd->nbonded_global += nexcl;
672 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0 && fplog)
673 {
674 fprintf(fplog,"There are %d inter charge-group exclusions,\n"
675 "will use an extra communication step for exclusion forces for %s\n",
676 dd->n_intercg_excl,eel_names[ir->coulombtype]);
677 }
678 }
679
680 natoms = mtop->natoms;
681
682 if (vsite && vsite->n_intercg_vsite > 0)
683 {
684 if (fplog)
685 {
686 fprintf(fplog,"There are %d inter charge-group virtual sites,\n"
687 "will an extra communication step for selected coordinates and forces\n",
688 vsite->n_intercg_vsite);
689 }
690 init_domdec_vsites(dd,natoms);
691 }
692
693 if (dd->bInterCGcons)
694 {
695 init_domdec_constraints(dd,natoms,mtop,constr);
696 }
697 if (fplog)
698 {
699 fprintf(fplog,"\n");
700 }
701 }
702
703 static inline void add_ifunc(int nral,t_iatom *tiatoms,t_ilist *il)
704 {
705 t_iatom *liatoms;
706 int k;
707
708 if (il->nr+1+nral > il->nalloc)
709 {
710 il->nalloc += over_alloc_large(il->nr+1+nral);
711 srenew(il->iatoms,il->nalloc);
712 }
713 liatoms = il->iatoms + il->nr;
714 for(k=0; k<=nral; k++)
715 {
716 liatoms[k] = tiatoms[k];
717 }
718 il->nr += 1 + nral;
719 }
720
721 static void add_posres(int mol,int a_mol,gmx_molblock_t *molb,
722 t_iatom *iatoms,t_idef *idef)
723 {
724 int n,a_molb;
725 t_iparams *ip;
726
727 /* This position restraint has not been added yet,
728 * so it's index is the current number of position restraints.
729 */
730 n = idef->il[F_POSRES].nr/2;
731 if (n+1 > idef->iparams_posres_nalloc)
732 {
733 idef->iparams_posres_nalloc = over_alloc_dd(n+1);
734 srenew(idef->iparams_posres,idef->iparams_posres_nalloc);
735 }
736 ip = &idef->iparams_posres[n];
737 /* Copy the force constants */
738 *ip = idef->iparams[iatoms[0]];
739
740 /* Get the position restriant coordinats from the molblock */
741 a_molb = mol*molb->natoms_mol + a_mol;
742 if (a_molb >= molb->nposres_xA)
743 {
744 gmx_incons("Not enough position restraint coordinates");
745 }
746 ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
747 ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
748 ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
749 if (molb->nposres_xB > 0)
750 {
751 ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
752 ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
753 ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
754 }
755 else
756 {
757 ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
758 ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
759 ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
760 }
761 /* Set the parameter index for idef->iparams_posre */
762 iatoms[0] = n;
763 }
764
765 static void add_vsite(gmx_ga2la_t ga2la,int *index,int *rtil,
766 int ftype,int nral,
767 bool bHomeA,int a,int a_gl,int a_mol,
768 t_iatom *iatoms,
769 t_idef *idef,int **vsite_pbc,int *vsite_pbc_nalloc)
770 {
771 int k,ak_gl,vsi,pbc_a_mol;
772 t_iatom tiatoms[1+MAXATOMLIST],*iatoms_r;
773 int j,ftype_r,nral_r;
774
775 /* Copy the type */
776 tiatoms[0] = iatoms[0];
777
778 if (bHomeA)
779 {
780 /* We know the local index of the first atom */
781 tiatoms[1] = a;
782 }
783 else
784 {
785 /* Convert later in make_local_vsites */
786 tiatoms[1] = -a_gl - 1;
787 }
788
789 for(k=2; k<1+nral; k++)
790 {
791 ak_gl = a_gl + iatoms[k] - a_mol;
792 if (!ga2la_get_home(ga2la,ak_gl,&tiatoms[k]))
793 {
794 /* Copy the global index, convert later in make_local_vsites */
795 tiatoms[k] = -(ak_gl + 1);
796 }
797 }
798
799 /* Add this interaction to the local topology */
800 add_ifunc(nral,tiatoms,&idef->il[ftype]);
801 if (vsite_pbc)
802 {
803 vsi = idef->il[ftype].nr/(1+nral) - 1;
804 if (vsi >= vsite_pbc_nalloc[ftype-F_VSITE2])
805 {
806 vsite_pbc_nalloc[ftype-F_VSITE2] = over_alloc_large(vsi+1);
807 srenew(vsite_pbc[ftype-F_VSITE2],vsite_pbc_nalloc[ftype-F_VSITE2]);
808 }
809 if (bHomeA)
810 {
811 pbc_a_mol = iatoms[1+nral+1];
812 if (pbc_a_mol < 0)
813 {
814 /* The pbc flag is one of the following two options:
815 * -2: vsite and all constructing atoms are within the same cg, no pbc
816 * -1: vsite and its first constructing atom are in the same cg, do pbc
817 */
818 vsite_pbc[ftype-F_VSITE2][vsi] = pbc_a_mol;
819 }
820 else
821 {
822 /* Set the pbc atom for this vsite so we can make its pbc
823 * identical to the rest of the atoms in its charge group.
824 * Since the order of the atoms does not change within a charge
825 * group, we do not need the global to local atom index.
826 */
827 vsite_pbc[ftype-F_VSITE2][vsi] = a + pbc_a_mol - iatoms[1];
828 }
829 }
830 else
831 {
832 /* This vsite is non-home (required for recursion),
833 * and therefore there is no charge group to match pbc with.
834 * But we always turn on full_pbc to assure that higher order
835 * recursion works correctly.
836 */
837 vsite_pbc[ftype-F_VSITE2][vsi] = -1;
838 }
839 }
840
841 if (iatoms[1+nral])
842 {
843 /* Check for recursion */
844 for(k=2; k<1+nral; k++)
845 {
846 if ((iatoms[1+nral] & (2<<k)) && (tiatoms[k] < 0))
847 {
848 /* This construction atoms is a vsite and not a home atom */
849 if (gmx_debug_at)
850 {
851 fprintf(debug,"Constructing atom %d of vsite atom %d is a vsite and non-home\n",iatoms[k]+1,a_mol+1);
852 }
853 /* Find the vsite construction */
854
855 /* Check all interactions assigned to this atom */
856 j = index[iatoms[k]];
857 while (j < index[iatoms[k]+1])
858 {
859 ftype_r = rtil[j++];
860 nral_r = NRAL(ftype_r);
861 if (interaction_function[ftype_r].flags & IF_VSITE)
862 {
863 /* Add this vsite (recursion) */
864 add_vsite(ga2la,index,rtil,ftype_r,nral_r,
865 FALSE,-1,a_gl+iatoms[k]-iatoms[1],iatoms[k],
866 rtil+j,idef,vsite_pbc,vsite_pbc_nalloc);
867 j += 1 + nral_r + 2;
868 }
869 else
870 {
871 j += 1 + nral_r;
872 }
873 }
874 }
875 }
876 }
877 }
878
879 static void make_la2lc(gmx_domdec_t *dd)
880 {
881 int *cgindex,*la2lc,cg,a;
882
883 cgindex = dd->cgindex;
884
885 if (dd->nat_tot > dd->la2lc_nalloc)
886 {
887 dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
888 snew(dd->la2lc,dd->la2lc_nalloc);
889 }
890 la2lc = dd->la2lc;
891
892 /* Make the local atom to local cg index */
893 for(cg=0; cg<dd->ncg_tot; cg++)
894 {
895 for(a=cgindex[cg]; a<cgindex[cg+1]; a++)
896 {
897 la2lc[a] = cg;
898 }
899 }
900 }
901
902 static real dd_dist2(t_pbc *pbc_null,rvec *cg_cm,const int *la2lc,int i,int j)
903 {
904 rvec dx;
905
906 if (pbc_null)
907 {
908 pbc_dx_aiuc(pbc_null,cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
909 }
910 else
911 {
912 rvec_sub(cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
913 }
914
915 return norm2(dx);
916 }
917
918 static int make_local_bondeds(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
919 gmx_molblock_t *molb,
920 bool bRCheckMB,ivec rcheck,bool bRCheck2B,
921 real rc,
922 int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
923 t_idef *idef,gmx_vsite_t *vsite)
924 {
925 int nzone,nizone,ic,la0,la1,i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,d,k;
926 int *index,*rtil,**vsite_pbc,*vsite_pbc_nalloc;
927 t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
928 bool bBCheck,bUse,bLocal;
929 real rc2;
930 ivec k_zero,k_plus;
931 gmx_ga2la_t ga2la;
932 int a_loc;
933 int kc;
934 gmx_domdec_ns_ranges_t *izone;
935 gmx_reverse_top_t *rt;
936 gmx_molblock_ind_t *mbi;
937 int nbonded_local;
938
939 nzone = zones->n;
940 nizone = zones->nizone;
941 izone = zones->izone;
942
943 rc2 = rc*rc;
944
945 if (vsite && vsite->n_intercg_vsite > 0)
946 {
947 vsite_pbc = vsite->vsite_pbc_loc;
948 vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
949 }
950 else
951 {
952 vsite_pbc = NULL;
953 vsite_pbc_nalloc = NULL;
954 }
955
956 rt = dd->reverse_top;
957
958 bBCheck = rt->bBCheck;
959
960 /* Clear the counts */
961 for(ftype=0; ftype<F_NRE; ftype++)
962 {
963 idef->il[ftype].nr = 0;
964 }
965 nbonded_local = 0;
966
967 mbi = rt->mbi;
968
969 ga2la = dd->ga2la;
970
971 for(ic=0; ic<nzone; ic++)
972 {
973 la0 = dd->cgindex[zones->cg_range[ic]];
974 la1 = dd->cgindex[zones->cg_range[ic+1]];
975 for(i=la0; i<la1; i++)
976 {
977 /* Get the global atom number */
978 i_gl = dd->gatindex[i];
979 global_atomnr_to_moltype_ind(mbi,i_gl,&mb,&mt,&mol,&i_mol);
980 /* Check all interactions assigned to this atom */
981 index = rt->ril_mt[mt].index;
982 rtil = rt->ril_mt[mt].il;
983 j = index[i_mol];
984 while (j < index[i_mol+1])
985 {
986 ftype = rtil[j++];
987 iatoms = rtil + j;
988 nral = NRAL(ftype);
989 if (interaction_function[ftype].flags & IF_VSITE)
990 {
991 /* The vsite construction goes where the vsite itself is */
992 if (ic == 0)
993 {
994 add_vsite(dd->ga2la,index,rtil,ftype,nral,
995 TRUE,i,i_gl,i_mol,
996 iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
997 }
998 j += 1 + nral + 2;
999 }
1000 else
1001 {
1002 /* Copy the type */
1003 tiatoms[0] = iatoms[0];
1004
1005 if (nral == 1)
1006 {
1007 /* Assign single-body interactions to the home zone */
1008 if (ic == 0)
1009 {
1010 bUse = TRUE;
1011 tiatoms[1] = i;
1012 if (ftype == F_POSRES)
1013 {
1014 add_posres(mol,i_mol,&molb[mb],tiatoms,idef);
1015 }
1016 }
1017 else
1018 {
1019 bUse = FALSE;
1020 }
1021 }
1022 else if (nral == 2)
1023 {
1024 /* This is a two-body interaction, we can assign
1025 * analogous to the non-bonded assignments.
1026 */
1027 if (!ga2la_get(ga2la,i_gl+iatoms[2]-i_mol,&a_loc,&kc))
1028 {
1029 bUse = FALSE;
1030 }
1031 else
1032 {
1033 if (kc >= nzone)
1034 {
1035 kc -= nzone;
1036 }
1037 /* Check zone interaction assignments */
1038 bUse = ((ic < nizone && ic <= kc &&
1039 izone[ic].j0 <= kc && kc < izone[ic].j1) ||
1040 (kc < nizone && ic > kc &&
1041 izone[kc].j0 <= ic && ic < izone[kc].j1));
1042 if (bUse)
1043 {
1044 tiatoms[1] = i;
1045 tiatoms[2] = a_loc;
1046 /* If necessary check the cgcm distance */
1047 if (bRCheck2B &&
1048 dd_dist2(pbc_null,cg_cm,la2lc,
1049 tiatoms[1],tiatoms[2]) >= rc2)
1050 {
1051 bUse = FALSE;
1052 }
1053 }
1054 }
1055 }
1056 else
1057 {
1058 /* Assign this multi-body bonded interaction to
1059 * the local node if we have all the atoms involved
1060 * (local or communicated) and the minimum zone shift
1061 * in each dimension is zero, for dimensions
1062 * with 2 DD cells an extra check may be necessary.
1063 */
1064 bUse = TRUE;
1065 clear_ivec(k_zero);
1066 clear_ivec(k_plus);
1067 for(k=1; k<=nral && bUse; k++)
1068 {
1069 bLocal = ga2la_get(ga2la,i_gl+iatoms[k]-i_mol,
1070 &a_loc,&kc);
1071 if (!bLocal || kc >= zones->n)
1072 {
1073 /* We do not have this atom of this interaction
1074 * locally, or it comes from more than one cell
1075 * away.
1076 */
1077 bUse = FALSE;
1078 }
1079 else
1080 {
1081 tiatoms[k] = a_loc;
1082 for(d=0; d<DIM; d++)
1083 {
1084 if (zones->shift[kc][d] == 0)
1085 {
1086 k_zero[d] = k;
1087 }
1088 else
1089 {
1090 k_plus[d] = k;
1091 }
1092 }
1093 }
1094 }
1095 bUse = (bUse &&
1096 k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
1097 if (bRCheckMB)
1098 {
1099 for(d=0; (d<DIM && bUse); d++)
1100 {
1101 /* Check if the cg_cm distance falls within
1102 * the cut-off to avoid possible multiple
1103 * assignments of bonded interactions.
1104 */
1105 if (rcheck[d] &&
1106 k_plus[d] &&
1107 dd_dist2(pbc_null,cg_cm,la2lc,
1108 tiatoms[k_zero[d]],tiatoms[k_plus[d]]) >= rc2)
1109 {
1110 bUse = FALSE;
1111 }
1112 }
1113 }
1114 }
1115 if (bUse)
1116 {
1117 /* Add this interaction to the local topology */
1118 add_ifunc(nral,tiatoms,&idef->il[ftype]);
1119 /* Sum so we can check in global_stat
1120 * if we have everything.
1121 */
1122 if (bBCheck ||
1123 !(interaction_function[ftype].flags & IF_LIMZERO))
1124 {
1125 nbonded_local++;
1126 }
1127 }
1128 j += 1 + nral;
1129 }
1130 }
1131 }
1132 }
1133
1134 return nbonded_local;
1135 }
1136
1137 static int make_local_bondeds_intracg(gmx_domdec_t *dd,gmx_molblock_t *molb,
1138 t_idef *idef,gmx_vsite_t *vsite)
1139 {
1140 int i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,k;
1141 int *index,*rtil,**vsite_pbc,*vsite_pbc_nalloc;
1142 t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
1143 gmx_reverse_top_t *rt;
1144 gmx_molblock_ind_t *mbi;
1145 int nbonded_local;
1146
1147 if (vsite && vsite->n_intercg_vsite > 0)
1148 {
1149 vsite_pbc = vsite->vsite_pbc_loc;
1150 vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
1151 }
1152 else
1153 {
1154 vsite_pbc = NULL;
1155 vsite_pbc_nalloc = NULL;
1156 }
1157
1158 /* Clear the counts */
1159 for(ftype=0; ftype<F_NRE; ftype++)
1160 {
1161 idef->il[ftype].nr = 0;
1162 }
1163 nbonded_local = 0;
1164
1165 rt = dd->reverse_top;
1166
1167 if (rt->ril_mt_tot_size == 0)
1168 {
1169 /* There are no interactions to assign */
1170 return nbonded_local;
1171 }
1172
1173 mbi = rt->mbi;
1174
1175 for(i=0; i<dd->nat_home; i++)
1176 {
1177 /* Get the global atom number */
1178 i_gl = dd->gatindex[i];
1179 global_atomnr_to_moltype_ind(mbi,i_gl,&mb,&mt,&mol,&i_mol);
1180 /* Check all interactions assigned to this atom */
1181 index = rt->ril_mt[mt].index;
1182 rtil = rt->ril_mt[mt].il;
1183 /* Check all interactions assigned to this atom */
1184 j = index[i_mol];
1185 while (j < index[i_mol+1])
1186 {
1187 ftype = rtil[j++];
1188 iatoms = rtil + j;
1189 nral = NRAL(ftype);
1190 if (interaction_function[ftype].flags & IF_VSITE)
1191 {
1192 /* The vsite construction goes where the vsite itself is */
1193 add_vsite(dd->ga2la,index,rtil,ftype,nral,
1194 TRUE,i,i_gl,i_mol,
1195 iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
1196 j += 1 + nral + 2;
1197 }
1198 else
1199 {
1200 /* Copy the type */
1201 tiatoms[0] = iatoms[0];
1202 tiatoms[1] = i;
1203 for(k=2; k<=nral; k++)
1204 {
1205 tiatoms[k] = i + iatoms[k] - iatoms[1];
1206 }
1207 if (ftype == F_POSRES)
1208 {
1209 add_posres(mol,i_mol,&molb[mb],tiatoms,idef);
1210 }
1211 /* Add this interaction to the local topology */
1212 add_ifunc(nral,tiatoms,&idef->il[ftype]);
1213 /* Sum so we can check in global_stat if we have everything */
1214 nbonded_local++;
1215 j += 1 + nral;
1216 }
1217 }
1218 }
1219
1220 return nbonded_local;
1221 }
1222
1223 static int make_local_exclusions(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1224 gmx_mtop_t *mtop,
1225 bool bRCheck,real rc,
1226 int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1227 t_forcerec *fr,
1228 t_blocka *lexcls)
1229 {
1230 int nizone,n,count,ic,jla0,jla1,jla;
1231 int cg,la0,la1,la,a_gl,mb,mt,mol,a_mol,j,aj_mol;
1232 t_blocka *excls;
1233 gmx_ga2la_t ga2la;
1234 int a_loc;
1235 int cell;
1236 gmx_molblock_ind_t *mbi;
1237 real rc2;
1238
1239 /* Since for RF and PME we need to loop over the exclusions
1240 * we should store each exclusion only once. This is done
1241 * using the same zone scheme as used for neighbor searching.
1242 * The exclusions involving non-home atoms are stored only
1243 * one way: atom j is in the excl list of i only for j > i,
1244 * where i and j are local atom numbers.
1245 */
1246
1247 lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1248 if (lexcls->nr+1 > lexcls->nalloc_index)
1249 {
1250 lexcls->nalloc_index = over_alloc_dd(lexcls->nr)+1;
1251 srenew(lexcls->index,lexcls->nalloc_index);
1252 }
1253
1254 mbi = dd->reverse_top->mbi;
1255
1256 ga2la = dd->ga2la;
1257
1258 rc2 = rc*rc;
1259
1260 if (dd->n_intercg_excl)
1261 {
1262 nizone = zones->nizone;
1263 }
1264 else
1265 {
1266 nizone = 1;
1267 }
1268 n = 0;
1269 count = 0;
1270 for(ic=0; ic<nizone; ic++)
1271 {
1272 jla0 = dd->cgindex[zones->izone[ic].jcg0];
1273 jla1 = dd->cgindex[zones->izone[ic].jcg1];
1274 for(cg=zones->cg_range[ic]; cg<zones->cg_range[ic+1]; cg++)
1275 {
1276 /* Here we assume the number of exclusions in one charge group
1277 * is never larger than 1000.
1278 */
1279 if (n+1000 > lexcls->nalloc_a)
1280 {
1281 lexcls->nalloc_a = over_alloc_large(n+1000);
1282 srenew(lexcls->a,lexcls->nalloc_a);
1283 }
1284 la0 = dd->cgindex[cg];
1285 la1 = dd->cgindex[cg+1];
1286 if (GET_CGINFO_EXCL_INTER(fr->cginfo[cg]) ||
1287 !GET_CGINFO_EXCL_INTRA(fr->cginfo[cg]))
1288 {
1289 /* Copy the exclusions from the global top */
1290 for(la=la0; la<la1; la++) {
1291 lexcls->index[la] = n;
1292 a_gl = dd->gatindex[la];
1293 global_atomnr_to_moltype_ind(mbi,a_gl,&mb,&mt,&mol,&a_mol);
1294 excls = &mtop->moltype[mt].excls;
1295 for(j=excls->index[a_mol]; j<excls->index[a_mol+1]; j++)
1296 {
1297 aj_mol = excls->a[j];
1298 /* This computation of jla is only correct intra-cg */
1299 jla = la + aj_mol - a_mol;
1300 if (jla >= la0 && jla < la1)
1301 {
1302 /* This is an intra-cg exclusion. We can skip
1303 * the global indexing and distance checking.
1304 */
1305 /* Intra-cg exclusions are only required
1306 * for the home zone.
1307 */
1308 if (ic == 0)
1309 {
1310 lexcls->a[n++] = jla;
1311 /* Check to avoid double counts */
1312 if (jla > la)
1313 {
1314 count++;
1315 }
1316 }
1317 }
1318 else
1319 {
1320 /* This is a inter-cg exclusion */
1321 /* Since exclusions are pair interactions,
1322 * just like non-bonded interactions,
1323 * they can be assigned properly up
1324 * to the DD cutoff (not cutoff_min as
1325 * for the other bonded interactions).
1326 */
1327 if (ga2la_get(ga2la,a_gl+aj_mol-a_mol,&jla,&cell))
1328 {
1329 if (ic == 0 && cell == 0)
1330 {
1331 lexcls->a[n++] = jla;
1332 /* Check to avoid double counts */
1333 if (jla > la)
1334 {
1335 count++;
1336 }
1337 }
1338 else if (jla >= jla0 && jla < jla1 &&
1339 (!bRCheck ||
1340 dd_dist2(pbc_null,cg_cm,la2lc,la,jla) < rc2))
1341 {
1342 /* jla > la, since jla0 > la */
1343 lexcls->a[n++] = jla;
1344 count++;
1345 }
1346 }
1347 }
1348 }
1349 }
1350 }
1351 else
1352 {
1353 /* There are no inter-cg excls and this cg is self-excluded.
1354 * These exclusions are only required for zone 0,
1355 * since other zones do not see themselves.
1356 */
1357 if (ic == 0)
1358 {
1359 for(la=la0; la<la1; la++)
1360 {
1361 lexcls->index[la] = n;
1362 for(j=la0; j<la1; j++)
1363 {
1364 lexcls->a[n++] = j;
1365 }
1366 }
1367 count += ((la1 - la0)*(la1 - la0 - 1))/2;
1368 }
1369 else
1370 {
1371 /* We don't need exclusions for this cg */
1372 for(la=la0; la<la1; la++)
1373 {
1374 lexcls->index[la] = n;
1375 }
1376 }
1377 }
1378 }
1379 }
1380 if (dd->n_intercg_excl == 0)
1381 {
1382 /* There are no exclusions involving non-home charge groups,
1383 * but we need to set the indices for neighborsearching.
1384 */
1385 la0 = dd->cgindex[zones->izone[0].cg1];
1386 for(la=la0; la<lexcls->nr; la++)
1387 {
1388 lexcls->index[la] = n;
1389 }
1390 }
1391 lexcls->index[lexcls->nr] = n;
1392 lexcls->nra = n;
1393 if (dd->n_intercg_excl == 0)
1394 {
1395 /* nr is only used to loop over the exclusions for Ewald and RF,
1396 * so we can set it to the number of home atoms for efficiency.
1397 */
1398 lexcls->nr = dd->cgindex[zones->izone[0].cg1];
1399 }
1400 if (debug)
1401 {
1402 fprintf(debug,"We have %d exclusions, check count %d\n",
1403 lexcls->nra,count);
1404 }
1405
1406 return count;
1407 }
1408
1409 void dd_make_local_cgs(gmx_domdec_t *dd,t_block *lcgs)
1410 {
1411 lcgs->nr = dd->ncg_tot;
1412 lcgs->index = dd->cgindex;
1413 }
1414
1415 void dd_make_local_top(FILE *fplog,
1416 gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1417 int npbcdim,matrix box,
1418 rvec cellsize_min,ivec npulse,
1419 t_forcerec *fr,gmx_vsite_t *vsite,
1420 gmx_mtop_t *mtop,gmx_localtop_t *ltop)
1421 {
1422 bool bUniqueExcl,bRCheckMB,bRCheck2B,bRCheckExcl;
1423 real rc=-1;
1424 ivec rcheck;
1425 int d,nexcl;
1426 t_pbc pbc,*pbc_null=NULL;
1427
1428 if (debug)
1429 {
1430 fprintf(debug,"Making local topology\n");
1431 }
1432
1433 dd_make_local_cgs(dd,&ltop->cgs);
1434
1435 bRCheckMB = FALSE;
1436 bRCheck2B = FALSE;
1437 bRCheckExcl = FALSE;
1438
1439 if (!dd->reverse_top->bMultiCGmols)
1440 {
1441 /* We don't need checks, assign all interactions with local atoms */
1442
1443 dd->nbonded_local = make_local_bondeds_intracg(dd,mtop->molblock,
1444 &ltop->idef,vsite);
1445 }
1446 else
1447 {
1448 /* We need to check to which cell bondeds should be assigned */
1449 rc = dd_cutoff_twobody(dd);
1450 if (debug)
1451 {
1452 fprintf(debug,"Two-body bonded cut-off distance is %g\n",rc);
1453 }
1454
1455 /* Should we check cg_cm distances when assigning bonded interactions? */
1456 for(d=0; d<DIM; d++)
1457 {
1458 rcheck[d] = FALSE;
1459 /* Only need to check for dimensions where the part of the box
1460 * that is not communicated is smaller than the cut-off.
1461 */
1462 if (d < npbcdim && dd->nc[d] > 1 &&
1463 (dd->nc[d] - npulse[d])*cellsize_min[d] < 2*rc)
1464 {
1465 if (dd->nc[d] == 2)
1466 {
1467 rcheck[d] = TRUE;
1468 bRCheckMB = TRUE;
1469 }
1470 /* Check for interactions between two atoms,
1471 * where we can allow interactions up to the cut-off,
1472 * instead of up to the smallest cell dimension.
1473 */
1474 bRCheck2B = TRUE;
1475 }
1476 if (debug)
1477 {
1478 fprintf(debug,
1479 "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
1480 d,cellsize_min[d],d,rcheck[d],bRCheck2B);
1481 }
1482 }
1483 if (dd->reverse_top->bExclRequired)
1484 {
1485 bRCheckExcl = bRCheck2B;
1486 }
1487 else
1488 {
1489 /* If we don't have forces on exclusions,
1490 * we don't care about exclusions being assigned mulitple times.
1491 */
1492 bRCheckExcl = FALSE;
1493 }
1494 if (bRCheckMB || bRCheck2B)
1495 {
1496 make_la2lc(dd);
1497 if (fr->bMolPBC)
1498 {
1499 set_pbc_dd(&pbc,fr->ePBC,dd,TRUE,box);
1500 pbc_null = &pbc;
1501 }
1502 else
1503 {
1504 pbc_null = NULL;
1505 }
1506 }
1507
1508 dd->nbonded_local = make_local_bondeds(dd,zones,mtop->molblock,
1509 bRCheckMB,rcheck,bRCheck2B,rc,
1510 dd->la2lc,
1511 pbc_null,fr->cg_cm,
1512 &ltop->idef,vsite);
1513 }
1514
1515 /* The ilist is not sorted yet,
1516 * we can only do this when we have the charge arrays.
1517 */
1518 ltop->idef.ilsort = ilsortUNKNOWN;
1519
1520 nexcl = make_local_exclusions(dd,zones,mtop,bRCheckExcl,
1521 rc,dd->la2lc,pbc_null,fr->cg_cm,
1522 fr,&ltop->excls);
1523
1524 if (dd->reverse_top->bExclRequired)
1525 {
1526 dd->nbonded_local += nexcl;
1527 }
1528
1529 ltop->atomtypes = mtop->atomtypes;
1530
1531 /* For an error message only */
1532 dd->reverse_top->err_top_global = mtop;
1533 dd->reverse_top->err_top_local = ltop;
1534 }
1535
1536 void dd_sort_local_top(gmx_domdec_t *dd,t_mdatoms *mdatoms,
1537 gmx_localtop_t *ltop)
1538 {
1539 if (dd->reverse_top->ilsort == ilsortNO_FE)
1540 {
1541 ltop->idef.ilsort = ilsortNO_FE;
1542 }
1543 else
1544 {
1545 gmx_sort_ilist_fe(&ltop->idef,mdatoms->chargeA,mdatoms->chargeB);
1546 }
1547 }
1548
1549 gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global)
1550 {
1551 gmx_localtop_t *top;
1552 int i;
1553
1554 snew(top,1);
1555
1556 top->idef.ntypes = top_global->ffparams.ntypes;
1557 top->idef.atnr = top_global->ffparams.atnr;
1558 top->idef.functype = top_global->ffparams.functype;
1559 top->idef.iparams = top_global->ffparams.iparams;
1560 top->idef.fudgeQQ = top_global->ffparams.fudgeQQ;
1561 top->idef.cmap_grid= top_global->ffparams.cmap_grid;
1562
1563 for(i=0; i<F_NRE; i++)
1564 {
1565 top->idef.il[i].iatoms = NULL;
1566 top->idef.il[i].nalloc = 0;
1567 }
1568 top->idef.ilsort = ilsortUNKNOWN;
1569
1570 return top;
1571 }
1572
1573 void dd_init_local_state(gmx_domdec_t *dd,
1574 t_state *state_global,t_state *state_local)
1575 {
1576 int i,j, buf[4];
1577
1578 if (DDMASTER(dd))
1579 {
1580 buf[0] = state_global->flags;
1581 buf[1] = state_global->ngtc;
1582 buf[2] = state_global->nnhpres;
1583 buf[3] = state_global->nhchainlength;
1584 }
1585 dd_bcast(dd,4*sizeof(int),buf);
1586
1587 init_state(state_local,0,buf[1],buf[2],buf[3]);
1588 state_local->flags = buf[0];
1589
1590 /* With Langevin Dynamics we need to make proper storage space
1591 * in the global and local state for the random numbers.
1592 */
1593 if (state_local->flags & (1<<estLD_RNG))
1594 {
1595 if (DDMASTER(dd) && state_global->nrngi > 1)
1596 {
1597 state_global->nrng = dd->nnodes*gmx_rng_n();
1598 srenew(state_global->ld_rng,state_global->nrng);
1599 }
1600 state_local->nrng = gmx_rng_n();
1601 snew(state_local->ld_rng,state_local->nrng);
1602 }
1603 if (state_local->flags & (1<<estLD_RNGI))
1604 {
1605 if (DDMASTER(dd) && state_global->nrngi > 1)
1606 {
1607 state_global->nrngi = dd->nnodes;
1608 srenew(state_global->ld_rngi,state_global->nrngi);
1609 }
1610 snew(state_local->ld_rngi,1);
1611 }
1612 }
1613
1614 static void check_link(t_blocka *link,int cg_gl,int cg_gl_j)
1615 {
1616 int k,aj;
1617 bool bFound;
1618
1619 bFound = FALSE;
1620 for(k=link->index[cg_gl]; k<link->index[cg_gl+1]; k++)
1621 {
1622 if (link->a[k] == cg_gl_j)
1623 {
1624 bFound = TRUE;
1625 }
1626 }
1627 if (!bFound)
1628 {
1629 /* Add this charge group link */
1630 if (link->index[cg_gl+1]+1 > link->nalloc_a)
1631 {
1632 link->nalloc_a = over_alloc_large(link->index[cg_gl+1]+1);
1633 srenew(link->a,link->nalloc_a);
1634 }
1635 link->a[link->index[cg_gl+1]] = cg_gl_j;
1636 link->index[cg_gl+1]++;
1637 }
1638 }
1639
1640 static int *make_at2cg(t_block *cgs)
1641 {
1642 int *at2cg,cg,a;
1643
1644 snew(at2cg,cgs->index[cgs->nr]);
1645 for(cg=0; cg<cgs->nr; cg++)
1646 {
1647 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1648 {
1649 at2cg[a] = cg;
1650 }
1651 }
1652
1653 return at2cg;
1654 }
1655
1656 t_blocka *make_charge_group_links(gmx_mtop_t *mtop,gmx_domdec_t *dd,
1657 cginfo_mb_t *cginfo_mb)
1658 {
1659 gmx_reverse_top_t *rt;
1660 int mb,cg_offset,cg,cg_gl,a,aj,i,j,ftype,nral,nlink_mol,mol,ncgi;
1661 gmx_molblock_t *molb;
1662 gmx_moltype_t *molt;
1663 t_block *cgs;
1664 t_blocka *excls;
1665 int *a2c;
1666 gmx_reverse_ilist_t ril;
1667 t_blocka *link;
1668 cginfo_mb_t *cgi_mb;
1669
1670 /* For each charge group make a list of other charge groups
1671 * in the system that a linked to it via bonded interactions
1672 * which are also stored in reverse_top.
1673 */
1674
1675 rt = dd->reverse_top;
1676
1677 snew(link,1);
1678 snew(link->index,ncg_mtop(mtop)+1);
1679 link->nalloc_a = 0;
1680 link->a = NULL;
1681
1682 link->index[0] = 0;
1683 cg_offset = 0;
1684 ncgi = 0;
1685 for(mb=0; mb<mtop->nmolblock; mb++)
1686 {
1687 molb = &mtop->molblock[mb];
1688 if (molb->nmol == 0)
1689 {
1690 continue;
1691 }
1692 molt = &mtop->moltype[molb->type];
1693 cgs = &molt->cgs;
1694 excls = &molt->excls;
1695 a2c = make_at2cg(cgs);
1696 /* Make a reverse ilist in which the interactions are linked
1697 * to all atoms, not only the first atom as in gmx_reverse_top.
1698 * The constraints are discarded here.
1699 */
1700 make_reverse_ilist(molt,NULL,FALSE,FALSE,TRUE,&ril);
1701
1702 cgi_mb = &cginfo_mb[mb];
1703
1704 for(cg=0; cg<cgs->nr; cg++)
1705 {
1706 cg_gl = cg_offset + cg;
1707 link->index[cg_gl+1] = link->index[cg_gl];
1708 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1709 {
1710 i = ril.index[a];
1711 while (i < ril.index[a+1])
1712 {
1713 ftype = ril.il[i++];
1714 nral = NRAL(ftype);
1715 /* Skip the ifunc index */
1716 i++;
1717 for(j=0; j<nral; j++)
1718 {
1719 aj = ril.il[i+j];
1720 if (a2c[aj] != cg)
1721 {
1722 check_link(link,cg_gl,cg_offset+a2c[aj]);
1723 }
1724 }
1725 i += nral_rt(ftype);
1726 }
1727 if (rt->bExclRequired)
1728 {
1729 /* Exclusions always go both ways */
1730 for(j=excls->index[a]; j<excls->index[a+1]; j++)
1731 {
1732 aj = excls->a[j];
1733 if (a2c[aj] != cg)
1734 {
1735 check_link(link,cg_gl,cg_offset+a2c[aj]);
1736 }
1737 }
1738 }
1739 }
1740 if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
1741 {
1742 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
1743 ncgi++;
1744 }
1745 }
1746 nlink_mol = link->index[cg_offset+cgs->nr] - link->index[cg_offset];
1747
1748 cg_offset += cgs->nr;
1749
1750 destroy_reverse_ilist(&ril);
1751 sfree(a2c);
1752
1753 if (debug)
1754 {
1755 fprintf(debug,"molecule type '%s' %d cgs has %d cg links through bonded interac.\n",*molt->name,cgs->nr,nlink_mol);
1756 }
1757
1758 if (molb->nmol > 1)
1759 {
1760 /* Copy the data for the rest of the molecules in this block */
1761 link->nalloc_a += (molb->nmol - 1)*nlink_mol;
1762 srenew(link->a,link->nalloc_a);
1763 for(mol=1; mol<molb->nmol; mol++)
1764 {
1765 for(cg=0; cg<cgs->nr; cg++)
1766 {
1767 cg_gl = cg_offset + cg;
1768 link->index[cg_gl+1] =
1769 link->index[cg_gl+1-cgs->nr] + nlink_mol;
1770 for(j=link->index[cg_gl]; j<link->index[cg_gl+1]; j++)
1771 {
1772 link->a[j] = link->a[j-nlink_mol] + cgs->nr;
1773 }
1774 if (link->index[cg_gl+1] - link->index[cg_gl] > 0 &&
1775 cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
1776 {
1777 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
1778 ncgi++;
1779 }
1780 }
1781 cg_offset += cgs->nr;
1782 }
1783 }
1784 }
1785
1786 if (debug)
1787 {
1788 fprintf(debug,"Of the %d charge groups %d are linked via bonded interactions\n",ncg_mtop(mtop),ncgi);
1789 }
1790
1791 return link;
1792 }
1793
1794 static void bonded_cg_distance_mol(gmx_moltype_t *molt,int *at2cg,
1795 bool bBCheck,bool bExcl,rvec *cg_cm,
1796 real *r_2b,int *ft2b,int *a2_1,int *a2_2,
1797 real *r_mb,int *ftmb,int *am_1,int *am_2)
1798 {
1799 int ftype,nral,i,j,ai,aj,cgi,cgj;
1800 t_ilist *il;
1801 t_blocka *excls;
1802 real r2_2b,r2_mb,rij2;
1803
1804 r2_2b = 0;
1805 r2_mb = 0;
1806 for(ftype=0; ftype<F_NRE; ftype++)
1807 {
1808 if (dd_check_ftype(ftype,bBCheck,FALSE))
1809 {
1810 il = &molt->ilist[ftype];
1811 nral = NRAL(ftype);
1812 if (nral > 1)
1813 {
1814 for(i=0; i<il->nr; i+=1+nral)
1815 {
1816 for(ai=0; ai<nral; ai++)
1817 {
1818 cgi = at2cg[il->iatoms[i+1+ai]];
1819 for(aj=0; aj<nral; aj++) {
1820 cgj = at2cg[il->iatoms[i+1+aj]];
1821 if (cgi != cgj)
1822 {
1823 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
1824 if (nral == 2 && rij2 > r2_2b)
1825 {
1826 r2_2b = rij2;
1827 *ft2b = ftype;
1828 *a2_1 = il->iatoms[i+1+ai];
1829 *a2_2 = il->iatoms[i+1+aj];
1830 }
1831 if (nral > 2 && rij2 > r2_mb)
1832 {
1833 r2_mb = rij2;
1834 *ftmb = ftype;
1835 *am_1 = il->iatoms[i+1+ai];
1836 *am_2 = il->iatoms[i+1+aj];
1837 }
1838 }
1839 }
1840 }
1841 }
1842 }
1843 }
1844 }
1845 if (bExcl)
1846 {
1847 excls = &molt->excls;
1848 for(ai=0; ai<excls->nr; ai++)
1849 {
1850 cgi = at2cg[ai];
1851 for(j=excls->index[ai]; j<excls->index[ai+1]; j++) {
1852 cgj = at2cg[excls->a[j]];
1853 if (cgi != cgj)
1854 {
1855 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
1856 if (rij2 > r2_2b)
1857 {
1858 r2_2b = rij2;
1859 }
1860 }
1861 }
1862 }
1863 }
1864
1865 *r_2b = sqrt(r2_2b);
1866 *r_mb = sqrt(r2_mb);
1867 }
1868
1869 static void get_cgcm_mol(gmx_moltype_t *molt,gmx_ffparams_t *ffparams,
1870 int ePBC,t_graph *graph,matrix box,
1871 gmx_vsite_t *vsite,
1872 rvec *x,rvec *xs,rvec *cg_cm)
1873 {
1874 int n,i;
1875
1876 if (ePBC != epbcNONE)
1877 {
1878 mk_mshift(NULL,graph,ePBC,box,x);
1879
1880 shift_x(graph,box,x,xs);
1881 /* By doing an extra mk_mshift the molecules that are broken
1882 * because they were e.g. imported from another software
1883 * will be made whole again. Such are the healing powers
1884 * of GROMACS.
1885 */
1886 mk_mshift(NULL,graph,ePBC,box,xs);
1887 }
1888 else
1889 {
1890 /* We copy the coordinates so the original coordinates remain
1891 * unchanged, just to be 100% sure that we do not affect
1892 * binary reproducability of simulations.
1893 */
1894 n = molt->cgs.index[molt->cgs.nr];
1895 for(i=0; i<n; i++)
1896 {
1897 copy_rvec(x[i],xs[i]);
1898 }
1899 }
1900
1901 if (vsite)
1902 {
1903 construct_vsites(NULL,vsite,xs,NULL,0.0,NULL,
1904 ffparams->iparams,molt->ilist,
1905 epbcNONE,TRUE,NULL,NULL,NULL);
1906 }
1907
1908 calc_cgcm(NULL,0,molt->cgs.nr,&molt->cgs,xs,cg_cm);
1909 }
1910
1911 static int have_vsite_molt(gmx_moltype_t *molt)
1912 {
1913 int i;
1914 bool bVSite;
1915
1916 bVSite = FALSE;
1917 for(i=0; i<F_NRE; i++)
1918 {
1919 if ((interaction_function[i].flags & IF_VSITE) &&
1920 molt->ilist[i].nr > 0) {
1921 bVSite = TRUE;
1922 }
1923 }
1924
1925 return bVSite;
1926 }
1927
1928 void dd_bonded_cg_distance(FILE *fplog,
1929 gmx_domdec_t *dd,gmx_mtop_t *mtop,
1930 t_inputrec *ir,rvec *x,matrix box,
1931 bool bBCheck,
1932 real *r_2b,real *r_mb)
1933 {
1934 bool bExclRequired;
1935 int mb,cg_offset,at_offset,*at2cg,mol;
1936 t_graph graph;
1937 gmx_vsite_t *vsite;
1938 gmx_molblock_t *molb;
1939 gmx_moltype_t *molt;
1940 rvec *xs,*cg_cm;
1941 real rmol_2b,rmol_mb;
1942 int ft2b=-1,a_2b_1=-1,a_2b_2=-1,ftmb=-1,a_mb_1=-1,a_mb_2=-1;
1943 int ftm2b=-1,amol_2b_1=-1,amol_2b_2=-1,ftmmb=-1,amol_mb_1=-1,amol_mb_2=-1;
1944
1945 bExclRequired = IR_EXCL_FORCES(*ir);
1946
1947 /* For gmx_vsite_t everything 0 should work (without pbc) */
1948 snew(vsite,1);
1949
1950 *r_2b = 0;
1951 *r_mb = 0;
1952 cg_offset = 0;
1953 at_offset = 0;
1954 for(mb=0; mb<mtop->nmolblock; mb++)
1955 {
1956 molb = &mtop->molblock[mb];
1957 molt = &mtop->moltype[molb->type];
1958 if (molt->cgs.nr == 1 || molb->nmol == 0)
1959 {
1960 cg_offset += molb->nmol*molt->cgs.nr;
1961 at_offset += molb->nmol*molt->atoms.nr;
1962 }
1963 else
1964 {
1965 if (ir->ePBC != epbcNONE)
1966 {
1967 mk_graph_ilist(NULL,molt->ilist,0,molt->atoms.nr,FALSE,FALSE,
1968 &graph);
1969 }
1970
1971 at2cg = make_at2cg(&molt->cgs);
1972 snew(xs,molt->atoms.nr);
1973 snew(cg_cm,molt->cgs.nr);
1974 for(mol=0; mol<molb->nmol; mol++)
1975 {
1976 get_cgcm_mol(molt,&mtop->ffparams,ir->ePBC,&graph,box,
1977 have_vsite_molt(molt) ? vsite : NULL,
1978 x+at_offset,xs,cg_cm);
1979
1980 bonded_cg_distance_mol(molt,at2cg,bBCheck,bExclRequired,cg_cm,
1981 &rmol_2b,&ftm2b,&amol_2b_1,&amol_2b_2,
1982 &rmol_mb,&ftmmb,&amol_mb_1,&amol_mb_2);
1983 if (rmol_2b > *r_2b)
1984 {
1985 *r_2b = rmol_2b;
1986 ft2b = ftm2b;
1987 a_2b_1 = at_offset + amol_2b_1;
1988 a_2b_2 = at_offset + amol_2b_2;
1989 }
1990 if (rmol_mb > *r_mb)
1991 {
1992 *r_mb = rmol_mb;
1993 ftmb = ftmmb;
1994 a_mb_1 = at_offset + amol_mb_1;
1995 a_mb_2 = at_offset + amol_mb_2;
1996 }
1997
1998 cg_offset += molt->cgs.nr;
1999 at_offset += molt->atoms.nr;
2000 }
2001 sfree(cg_cm);
2002 sfree(xs);
2003 sfree(at2cg);
2004 if (ir->ePBC != epbcNONE)
2005 {
2006 done_graph(&graph);
2007 }
2008 }
2009 }
2010
2011 sfree(vsite);
2012
2013 if (fplog && (ft2b >= 0 || ftmb >= 0))
2014 {
2015 fprintf(fplog,
2016 "Initial maximum inter charge-group distances:\n");
2017 if (ft2b >= 0)
2018 {
2019 fprintf(fplog,
2020 " two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2021 *r_2b,interaction_function[ft2b].longname,
2022 a_2b_1+1,a_2b_2+1);
2023 }
2024 if (ftmb >= 0)
2025 {
2026 fprintf(fplog,
2027 " multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2028 *r_mb,interaction_function[ftmb].longname,
2029 a_mb_1+1,a_mb_2+1);
2030 }
2031 }
2032 }
Atempt to write QMMM interface for Gamess-US