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