1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This file is part of Gromacs Copyright (c) 1991-2008
5 * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
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.
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
16 * Gnomes, ROck Monsters And Chili Sauce
27 #include "domdec_network.h"
32 #include "chargegroup.h"
33 #include "gmx_random.h"
35 #include "mtop_util.h"
38 #include "gmx_ga2la.h"
40 /* for dd_init_local_state */
41 #define NITEM_DD_INIT_LOCAL_STATE 5
44 int *index; /* Index for each atom into il */
45 int *il; /* ftype|type|a0|...|an|ftype|... */
46 } gmx_reverse_ilist_t;
55 typedef struct gmx_reverse_top {
56 gmx_bool bExclRequired; /* Do we require all exclusions to be assigned? */
57 gmx_bool bConstr; /* Are there constraints in this revserse top? */
58 gmx_bool bBCheck; /* All bonded interactions have to be assigned? */
59 gmx_bool bMultiCGmols; /* Are the multi charge-group molecules? */
60 gmx_reverse_ilist_t *ril_mt; /* Reverse ilist for all moltypes */
62 int ilsort; /* The sorting state of bondeds for free energy */
63 gmx_molblock_ind_t *mbi;
65 /* Pointers only used for an error message */
66 gmx_mtop_t *err_top_global;
67 gmx_localtop_t *err_top_local;
70 static int nral_rt(int ftype)
72 /* Returns the number of atom entries for il in gmx_reverse_top_t */
76 if (interaction_function[ftype].flags & IF_VSITE)
78 /* With vsites the reverse topology contains
79 * two extra entries for PBC.
87 static gmx_bool dd_check_ftype(int ftype,gmx_bool bBCheck,gmx_bool bConstr)
89 return (((interaction_function[ftype].flags & IF_BOND) &&
90 !(interaction_function[ftype].flags & IF_VSITE) &&
91 (bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO))) ||
93 (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)));
96 static void print_error_header(FILE *fplog,char *moltypename,int nprint)
98 fprintf(fplog, "\nMolecule type '%s'\n",moltypename);
99 fprintf(stderr,"\nMolecule type '%s'\n",moltypename);
101 "the first %d missing interactions, except for exclusions:\n",
104 "the first %d missing interactions, except for exclusions:\n",
108 static void print_missing_interactions_mb(FILE *fplog,t_commrec *cr,
109 gmx_reverse_top_t *rt,
111 gmx_reverse_ilist_t *ril,
112 int a_start,int a_end,
113 int nat_mol,int nmol,
116 int nril_mol,*assigned,*gatindex;
117 int ftype,ftype_j,nral,i,j_mol,j,k,a0,a0_mol,mol,a,a_gl;
123 nril_mol = ril->index[nat_mol];
124 snew(assigned,nmol*nril_mol);
126 gatindex = cr->dd->gatindex;
127 for(ftype=0; ftype<F_NRE; ftype++)
129 if (dd_check_ftype(ftype,rt->bBCheck,rt->bConstr))
132 il = &idef->il[ftype];
134 for(i=0; i<il->nr; i+=1+nral)
136 a0 = gatindex[ia[1]];
137 /* Check if this interaction is in
138 * the currently checked molblock.
140 if (a0 >= a_start && a0 < a_end)
142 mol = (a0 - a_start)/nat_mol;
143 a0_mol = (a0 - a_start) - mol*nat_mol;
144 j_mol = ril->index[a0_mol];
146 while (j_mol < ril->index[a0_mol+1] && !bFound)
148 j = mol*nril_mol + j_mol;
149 ftype_j = ril->il[j_mol];
150 /* Here we need to check if this interaction has
151 * not already been assigned, since we could have
152 * multiply defined interactions.
154 if (ftype == ftype_j && ia[0] == ril->il[j_mol+1] &&
157 /* Check the atoms */
159 for(a=0; a<nral; a++)
161 if (gatindex[ia[1+a]] !=
162 a_start + mol*nat_mol + ril->il[j_mol+2+a])
172 j_mol += 2 + nral_rt(ftype_j);
176 gmx_incons("Some interactions seem to be assigned multiple times");
184 gmx_sumi(nmol*nril_mol,assigned,cr);
188 for(mol=0; mol<nmol; mol++)
191 while (j_mol < nril_mol)
193 ftype = ril->il[j_mol];
195 j = mol*nril_mol + j_mol;
196 if (assigned[j] == 0 &&
197 !(interaction_function[ftype].flags & IF_VSITE))
199 if (DDMASTER(cr->dd))
203 print_error_header(fplog,moltypename,nprint);
205 fprintf(fplog, "%20s atoms",
206 interaction_function[ftype].longname);
207 fprintf(stderr,"%20s atoms",
208 interaction_function[ftype].longname);
209 for(a=0; a<nral; a++) {
210 fprintf(fplog, "%5d",ril->il[j_mol+2+a]+1);
211 fprintf(stderr,"%5d",ril->il[j_mol+2+a]+1);
219 fprintf(fplog, " global");
220 fprintf(stderr," global");
221 for(a=0; a<nral; a++)
223 fprintf(fplog, "%6d",
224 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
225 fprintf(stderr,"%6d",
226 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
228 fprintf(fplog, "\n");
229 fprintf(stderr,"\n");
237 j_mol += 2 + nral_rt(ftype);
244 static void print_missing_interactions_atoms(FILE *fplog,t_commrec *cr,
245 gmx_mtop_t *mtop,t_idef *idef)
247 int mb,a_start,a_end;
248 gmx_molblock_t *molb;
249 gmx_reverse_top_t *rt;
251 rt = cr->dd->reverse_top;
253 /* Print the atoms in the missing interactions per molblock */
255 for(mb=0; mb<mtop->nmolblock; mb++)
257 molb = &mtop->molblock[mb];
259 a_end = a_start + molb->nmol*molb->natoms_mol;
261 print_missing_interactions_mb(fplog,cr,rt,
262 *(mtop->moltype[molb->type].name),
263 &rt->ril_mt[molb->type],
264 a_start,a_end,molb->natoms_mol,
270 void dd_print_missing_interactions(FILE *fplog,t_commrec *cr,int local_count, gmx_mtop_t *top_global, t_state *state_local)
272 int ndiff_tot,cl[F_NRE],n,ndiff,rest_global,rest_local;
276 gmx_mtop_t *err_top_global;
277 gmx_localtop_t *err_top_local;
281 err_top_global = dd->reverse_top->err_top_global;
282 err_top_local = dd->reverse_top->err_top_local;
286 fprintf(fplog,"\nNot all bonded interactions have been properly assigned to the domain decomposition cells\n");
290 ndiff_tot = local_count - dd->nbonded_global;
292 for(ftype=0; ftype<F_NRE; ftype++)
295 cl[ftype] = err_top_local->idef.il[ftype].nr/(1+nral);
298 gmx_sumi(F_NRE,cl,cr);
302 fprintf(fplog,"\nA list of missing interactions:\n");
303 fprintf(stderr,"\nA list of missing interactions:\n");
304 rest_global = dd->nbonded_global;
305 rest_local = local_count;
306 for(ftype=0; ftype<F_NRE; ftype++)
308 /* In the reverse and local top all constraints are merged
309 * into F_CONSTR. So in the if statement we skip F_CONSTRNC
310 * and add these constraints when doing F_CONSTR.
312 if (((interaction_function[ftype].flags & IF_BOND) &&
313 (dd->reverse_top->bBCheck
314 || !(interaction_function[ftype].flags & IF_LIMZERO)))
316 || (dd->reverse_top->bConstr && ftype == F_CONSTR))
319 n = gmx_mtop_ftype_count(err_top_global,ftype);
320 if (ftype == F_CONSTR)
322 n += gmx_mtop_ftype_count(err_top_global,F_CONSTRNC);
324 ndiff = cl[ftype] - n;
327 sprintf(buf,"%20s of %6d missing %6d",
328 interaction_function[ftype].longname,n,-ndiff);
329 fprintf(fplog,"%s\n",buf);
330 fprintf(stderr,"%s\n",buf);
333 rest_local -= cl[ftype];
337 ndiff = rest_local - rest_global;
340 sprintf(buf,"%20s of %6d missing %6d","exclusions",
342 fprintf(fplog,"%s\n",buf);
343 fprintf(stderr,"%s\n",buf);
347 print_missing_interactions_atoms(fplog,cr,err_top_global,
348 &err_top_local->idef);
349 write_dd_pdb("dd_dump_err",0,"dump",top_global,cr,
350 -1,state_local->x,state_local->box);
355 gmx_incons("One or more interactions were multiple assigned in the domain decompostion");
359 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));
364 static void global_atomnr_to_moltype_ind(gmx_molblock_ind_t *mbi,int i_gl,
365 int *mb,int *mt,int *mol,int *i_mol)
370 while (i_gl >= mbi->a_end) {
376 *mol = (i_gl - mbi->a_start) / mbi->natoms_mol;
377 *i_mol = (i_gl - mbi->a_start) - (*mol)*mbi->natoms_mol;
380 static int count_excls(t_block *cgs,t_blocka *excls,int *n_intercg_excl)
382 int n,n_inter,cg,at0,at1,at,excl,atj;
386 for(cg=0; cg<cgs->nr; cg++)
388 at0 = cgs->index[cg];
389 at1 = cgs->index[cg+1];
390 for(at=at0; at<at1; at++)
392 for(excl=excls->index[at]; excl<excls->index[at+1]; excl++)
394 atj = excls->a[excl];
398 if (atj < at0 || atj >= at1)
410 static int low_make_reverse_ilist(t_ilist *il_mt,t_atom *atom,
413 gmx_bool bConstr,gmx_bool bBCheck,
414 int *r_index,int *r_il,
415 gmx_bool bLinkToAllAtoms,
418 int ftype,nral,i,j,nlink,link;
426 for(ftype=0; ftype<F_NRE; ftype++)
428 if ((interaction_function[ftype].flags & (IF_BOND | IF_VSITE)) ||
430 (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC))) {
431 bVSite = (interaction_function[ftype].flags & IF_VSITE);
435 for(i=0; i<il->nr; i+=1+nral)
442 /* We don't need the virtual sites for the cg-links */
452 /* Couple to the first atom in the interaction */
455 for(link=0; link<nlink; link++)
460 r_il[r_index[a]+count[a]] =
461 (ftype == F_CONSTRNC ? F_CONSTR : ftype);
462 r_il[r_index[a]+count[a]+1] = ia[0];
463 for(j=1; j<1+nral; j++)
465 /* Store the molecular atom number */
466 r_il[r_index[a]+count[a]+1+j] = ia[j];
469 if (interaction_function[ftype].flags & IF_VSITE)
473 /* Add an entry to iatoms for storing
474 * which of the constructing atoms are
477 r_il[r_index[a]+count[a]+2+nral] = 0;
478 for(j=2; j<1+nral; j++)
480 if (atom[ia[j]].ptype == eptVSite)
482 r_il[r_index[a]+count[a]+2+nral] |= (2<<j);
485 /* Store vsite pbc atom in a second extra entry */
486 r_il[r_index[a]+count[a]+2+nral+1] =
487 (vsite_pbc ? vsite_pbc[ftype-F_VSITE2][i/(1+nral)] : -2);
492 /* We do not count vsites since they are always
493 * uniquely assigned and can be assigned
494 * to multiple nodes with recursive vsites.
497 !(interaction_function[ftype].flags & IF_LIMZERO))
502 count[a] += 2 + nral_rt(ftype);
511 static int make_reverse_ilist(gmx_moltype_t *molt,
513 gmx_bool bConstr,gmx_bool bBCheck,
514 gmx_bool bLinkToAllAtoms,
515 gmx_reverse_ilist_t *ril_mt)
517 int nat_mt,*count,i,nint_mt;
519 /* Count the interactions */
520 nat_mt = molt->atoms.nr;
522 low_make_reverse_ilist(molt->ilist,molt->atoms.atom,vsite_pbc,
524 bConstr,bBCheck,NULL,NULL,
525 bLinkToAllAtoms,FALSE);
527 snew(ril_mt->index,nat_mt+1);
528 ril_mt->index[0] = 0;
529 for(i=0; i<nat_mt; i++)
531 ril_mt->index[i+1] = ril_mt->index[i] + count[i];
534 snew(ril_mt->il,ril_mt->index[nat_mt]);
536 /* Store the interactions */
538 low_make_reverse_ilist(molt->ilist,molt->atoms.atom,vsite_pbc,
541 ril_mt->index,ril_mt->il,
542 bLinkToAllAtoms,TRUE);
549 static void destroy_reverse_ilist(gmx_reverse_ilist_t *ril)
555 static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop,gmx_bool bFE,
556 int ***vsite_pbc_molt,
558 gmx_bool bBCheck,int *nint)
561 gmx_reverse_top_t *rt;
567 /* Should we include constraints (for SHAKE) in rt? */
568 rt->bConstr = bConstr;
569 rt->bBCheck = bBCheck;
571 rt->bMultiCGmols = FALSE;
572 snew(nint_mt,mtop->nmoltype);
573 snew(rt->ril_mt,mtop->nmoltype);
574 rt->ril_mt_tot_size = 0;
575 for(mt=0; mt<mtop->nmoltype; mt++)
577 molt = &mtop->moltype[mt];
578 if (molt->cgs.nr > 1)
580 rt->bMultiCGmols = TRUE;
583 /* Make the atom to interaction list for this molecule type */
585 make_reverse_ilist(molt,vsite_pbc_molt ? vsite_pbc_molt[mt] : NULL,
586 rt->bConstr,rt->bBCheck,FALSE,
589 rt->ril_mt_tot_size += rt->ril_mt[mt].index[molt->atoms.nr];
593 fprintf(debug,"The total size of the atom to interaction index is %d integers\n",rt->ril_mt_tot_size);
597 for(mb=0; mb<mtop->nmolblock; mb++)
599 *nint += mtop->molblock[mb].nmol*nint_mt[mtop->molblock[mb].type];
603 if (bFE && gmx_mtop_bondeds_free_energy(mtop))
605 rt->ilsort = ilsortFE_UNSORTED;
608 rt->ilsort = ilsortNO_FE;
611 /* Make a molblock index for fast searching */
612 snew(rt->mbi,mtop->nmolblock);
614 for(mb=0; mb<mtop->nmolblock; mb++)
616 rt->mbi[mb].a_start = i;
617 i += mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
618 rt->mbi[mb].a_end = i;
619 rt->mbi[mb].natoms_mol = mtop->molblock[mb].natoms_mol;
620 rt->mbi[mb].type = mtop->molblock[mb].type;
626 void dd_make_reverse_top(FILE *fplog,
627 gmx_domdec_t *dd,gmx_mtop_t *mtop,
628 gmx_vsite_t *vsite,gmx_constr_t constr,
629 t_inputrec *ir,gmx_bool bBCheck)
631 int mb,natoms,n_recursive_vsite,nexcl,nexcl_icg,a;
632 gmx_molblock_t *molb;
637 fprintf(fplog,"\nLinking all bonded interactions to atoms\n");
640 dd->reverse_top = make_reverse_top(mtop,ir->efep!=efepNO,
641 vsite ? vsite->vsite_pbc_molt : NULL,
643 bBCheck,&dd->nbonded_global);
645 if (dd->reverse_top->ril_mt_tot_size >= 200000 &&
647 mtop->nmolblock == 1 && mtop->molblock[0].nmol == 1)
649 /* mtop comes from a pre Gromacs 4 tpr file */
650 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";
653 fprintf(fplog,"\n%s\n\n",note);
657 fprintf(stderr,"\n%s\n\n",note);
661 dd->reverse_top->bExclRequired = IR_EXCL_FORCES(*ir);
664 dd->n_intercg_excl = 0;
665 for(mb=0; mb<mtop->nmolblock; mb++)
667 molb = &mtop->molblock[mb];
668 molt = &mtop->moltype[molb->type];
669 nexcl += molb->nmol*count_excls(&molt->cgs,&molt->excls,&nexcl_icg);
670 dd->n_intercg_excl += molb->nmol*nexcl_icg;
672 if (dd->reverse_top->bExclRequired)
674 dd->nbonded_global += nexcl;
675 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0 && fplog)
677 fprintf(fplog,"There are %d inter charge-group exclusions,\n"
678 "will use an extra communication step for exclusion forces for %s\n",
679 dd->n_intercg_excl,eel_names[ir->coulombtype]);
683 natoms = mtop->natoms;
685 if (vsite && vsite->n_intercg_vsite > 0)
689 fprintf(fplog,"There are %d inter charge-group virtual sites,\n"
690 "will an extra communication step for selected coordinates and forces\n",
691 vsite->n_intercg_vsite);
693 init_domdec_vsites(dd,natoms);
696 if (dd->bInterCGcons)
698 init_domdec_constraints(dd,natoms,mtop,constr);
706 static inline void add_ifunc(int nral,t_iatom *tiatoms,t_ilist *il)
711 if (il->nr+1+nral > il->nalloc)
713 il->nalloc += over_alloc_large(il->nr+1+nral);
714 srenew(il->iatoms,il->nalloc);
716 liatoms = il->iatoms + il->nr;
717 for(k=0; k<=nral; k++)
719 liatoms[k] = tiatoms[k];
724 static void add_posres(int mol,int a_mol,gmx_molblock_t *molb,
725 t_iatom *iatoms,t_idef *idef)
730 /* This position restraint has not been added yet,
731 * so it's index is the current number of position restraints.
733 n = idef->il[F_POSRES].nr/2;
734 if (n+1 > idef->iparams_posres_nalloc)
736 idef->iparams_posres_nalloc = over_alloc_dd(n+1);
737 srenew(idef->iparams_posres,idef->iparams_posres_nalloc);
739 ip = &idef->iparams_posres[n];
740 /* Copy the force constants */
741 *ip = idef->iparams[iatoms[0]];
743 /* Get the position restriant coordinats from the molblock */
744 a_molb = mol*molb->natoms_mol + a_mol;
745 if (a_molb >= molb->nposres_xA)
747 gmx_incons("Not enough position restraint coordinates");
749 ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
750 ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
751 ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
752 if (molb->nposres_xB > 0)
754 ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
755 ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
756 ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
760 ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
761 ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
762 ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
764 /* Set the parameter index for idef->iparams_posre */
768 static void add_fbposres(int mol,int a_mol,gmx_molblock_t *molb,
769 t_iatom *iatoms,t_idef *idef)
774 /* This flat-bottom position restraint has not been added yet,
775 * so it's index is the current number of position restraints.
777 n = idef->il[F_FBPOSRES].nr/2;
778 if (n+1 > idef->iparams_fbposres_nalloc)
780 idef->iparams_fbposres_nalloc = over_alloc_dd(n+1);
781 srenew(idef->iparams_fbposres,idef->iparams_fbposres_nalloc);
783 ip = &idef->iparams_fbposres[n];
784 /* Copy the force constants */
785 *ip = idef->iparams[iatoms[0]];
787 /* Get the position restriant coordinats from the molblock */
788 a_molb = mol*molb->natoms_mol + a_mol;
789 if (a_molb >= molb->nposres_xA)
791 gmx_incons("Not enough position restraint coordinates");
793 /* Take reference positions from A position of normal posres */
794 ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
795 ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
796 ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
798 /* Note: no B-type for flat-bottom posres */
800 /* Set the parameter index for idef->iparams_posre */
804 static void add_vsite(gmx_ga2la_t ga2la,int *index,int *rtil,
806 gmx_bool bHomeA,int a,int a_gl,int a_mol,
808 t_idef *idef,int **vsite_pbc,int *vsite_pbc_nalloc)
810 int k,ak_gl,vsi,pbc_a_mol;
811 t_iatom tiatoms[1+MAXATOMLIST],*iatoms_r;
812 int j,ftype_r,nral_r;
815 tiatoms[0] = iatoms[0];
819 /* We know the local index of the first atom */
824 /* Convert later in make_local_vsites */
825 tiatoms[1] = -a_gl - 1;
828 for(k=2; k<1+nral; k++)
830 ak_gl = a_gl + iatoms[k] - a_mol;
831 if (!ga2la_get_home(ga2la,ak_gl,&tiatoms[k]))
833 /* Copy the global index, convert later in make_local_vsites */
834 tiatoms[k] = -(ak_gl + 1);
838 /* Add this interaction to the local topology */
839 add_ifunc(nral,tiatoms,&idef->il[ftype]);
842 vsi = idef->il[ftype].nr/(1+nral) - 1;
843 if (vsi >= vsite_pbc_nalloc[ftype-F_VSITE2])
845 vsite_pbc_nalloc[ftype-F_VSITE2] = over_alloc_large(vsi+1);
846 srenew(vsite_pbc[ftype-F_VSITE2],vsite_pbc_nalloc[ftype-F_VSITE2]);
850 pbc_a_mol = iatoms[1+nral+1];
853 /* The pbc flag is one of the following two options:
854 * -2: vsite and all constructing atoms are within the same cg, no pbc
855 * -1: vsite and its first constructing atom are in the same cg, do pbc
857 vsite_pbc[ftype-F_VSITE2][vsi] = pbc_a_mol;
861 /* Set the pbc atom for this vsite so we can make its pbc
862 * identical to the rest of the atoms in its charge group.
863 * Since the order of the atoms does not change within a charge
864 * group, we do not need the global to local atom index.
866 vsite_pbc[ftype-F_VSITE2][vsi] = a + pbc_a_mol - iatoms[1];
871 /* This vsite is non-home (required for recursion),
872 * and therefore there is no charge group to match pbc with.
873 * But we always turn on full_pbc to assure that higher order
874 * recursion works correctly.
876 vsite_pbc[ftype-F_VSITE2][vsi] = -1;
882 /* Check for recursion */
883 for(k=2; k<1+nral; k++)
885 if ((iatoms[1+nral] & (2<<k)) && (tiatoms[k] < 0))
887 /* This construction atoms is a vsite and not a home atom */
890 fprintf(debug,"Constructing atom %d of vsite atom %d is a vsite and non-home\n",iatoms[k]+1,a_mol+1);
892 /* Find the vsite construction */
894 /* Check all interactions assigned to this atom */
895 j = index[iatoms[k]];
896 while (j < index[iatoms[k]+1])
899 nral_r = NRAL(ftype_r);
900 if (interaction_function[ftype_r].flags & IF_VSITE)
902 /* Add this vsite (recursion) */
903 add_vsite(ga2la,index,rtil,ftype_r,nral_r,
904 FALSE,-1,a_gl+iatoms[k]-iatoms[1],iatoms[k],
905 rtil+j,idef,vsite_pbc,vsite_pbc_nalloc);
918 static void make_la2lc(gmx_domdec_t *dd)
920 int *cgindex,*la2lc,cg,a;
922 cgindex = dd->cgindex;
924 if (dd->nat_tot > dd->la2lc_nalloc)
926 dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
927 snew(dd->la2lc,dd->la2lc_nalloc);
931 /* Make the local atom to local cg index */
932 for(cg=0; cg<dd->ncg_tot; cg++)
934 for(a=cgindex[cg]; a<cgindex[cg+1]; a++)
941 static real dd_dist2(t_pbc *pbc_null,rvec *cg_cm,const int *la2lc,int i,int j)
947 pbc_dx_aiuc(pbc_null,cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
951 rvec_sub(cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
957 static int make_local_bondeds(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
958 gmx_molblock_t *molb,
959 gmx_bool bRCheckMB,ivec rcheck,gmx_bool bRCheck2B,
961 int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
962 t_idef *idef,gmx_vsite_t *vsite)
964 int nzone,nizone,ic,la0,la1,i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,d,k;
965 int *index,*rtil,**vsite_pbc,*vsite_pbc_nalloc;
966 t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
967 gmx_bool bBCheck,bUse,bLocal;
973 gmx_domdec_ns_ranges_t *izone;
974 gmx_reverse_top_t *rt;
975 gmx_molblock_ind_t *mbi;
979 nizone = zones->nizone;
980 izone = zones->izone;
984 if (vsite && vsite->n_intercg_vsite > 0)
986 vsite_pbc = vsite->vsite_pbc_loc;
987 vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
992 vsite_pbc_nalloc = NULL;
995 rt = dd->reverse_top;
997 bBCheck = rt->bBCheck;
999 /* Clear the counts */
1000 for(ftype=0; ftype<F_NRE; ftype++)
1002 idef->il[ftype].nr = 0;
1010 for(ic=0; ic<nzone; ic++)
1012 la0 = dd->cgindex[zones->cg_range[ic]];
1013 la1 = dd->cgindex[zones->cg_range[ic+1]];
1014 for(i=la0; i<la1; i++)
1016 /* Get the global atom number */
1017 i_gl = dd->gatindex[i];
1018 global_atomnr_to_moltype_ind(mbi,i_gl,&mb,&mt,&mol,&i_mol);
1019 /* Check all interactions assigned to this atom */
1020 index = rt->ril_mt[mt].index;
1021 rtil = rt->ril_mt[mt].il;
1023 while (j < index[i_mol+1])
1028 if (interaction_function[ftype].flags & IF_VSITE)
1030 /* The vsite construction goes where the vsite itself is */
1033 add_vsite(dd->ga2la,index,rtil,ftype,nral,
1035 iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
1042 tiatoms[0] = iatoms[0];
1046 /* Assign single-body interactions to the home zone */
1051 if (ftype == F_POSRES)
1053 add_posres(mol,i_mol,&molb[mb],tiatoms,idef);
1055 else if (ftype == F_FBPOSRES)
1057 add_fbposres(mol,i_mol,&molb[mb],tiatoms,idef);
1067 /* This is a two-body interaction, we can assign
1068 * analogous to the non-bonded assignments.
1070 if (!ga2la_get(ga2la,i_gl+iatoms[2]-i_mol,&a_loc,&kc))
1080 /* Check zone interaction assignments */
1081 bUse = ((ic < nizone && ic <= kc &&
1082 izone[ic].j0 <= kc && kc < izone[ic].j1) ||
1083 (kc < nizone && ic > kc &&
1084 izone[kc].j0 <= ic && ic < izone[kc].j1));
1089 /* If necessary check the cgcm distance */
1091 dd_dist2(pbc_null,cg_cm,la2lc,
1092 tiatoms[1],tiatoms[2]) >= rc2)
1101 /* Assign this multi-body bonded interaction to
1102 * the local node if we have all the atoms involved
1103 * (local or communicated) and the minimum zone shift
1104 * in each dimension is zero, for dimensions
1105 * with 2 DD cells an extra check may be necessary.
1110 for(k=1; k<=nral && bUse; k++)
1112 bLocal = ga2la_get(ga2la,i_gl+iatoms[k]-i_mol,
1114 if (!bLocal || kc >= zones->n)
1116 /* We do not have this atom of this interaction
1117 * locally, or it comes from more than one cell
1125 for(d=0; d<DIM; d++)
1127 if (zones->shift[kc][d] == 0)
1139 k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
1142 for(d=0; (d<DIM && bUse); d++)
1144 /* Check if the cg_cm distance falls within
1145 * the cut-off to avoid possible multiple
1146 * assignments of bonded interactions.
1150 dd_dist2(pbc_null,cg_cm,la2lc,
1151 tiatoms[k_zero[d]],tiatoms[k_plus[d]]) >= rc2)
1160 /* Add this interaction to the local topology */
1161 add_ifunc(nral,tiatoms,&idef->il[ftype]);
1162 /* Sum so we can check in global_stat
1163 * if we have everything.
1166 !(interaction_function[ftype].flags & IF_LIMZERO))
1177 return nbonded_local;
1180 static int make_local_bondeds_intracg(gmx_domdec_t *dd,gmx_molblock_t *molb,
1181 t_idef *idef,gmx_vsite_t *vsite)
1183 int i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,k;
1184 int *index,*rtil,**vsite_pbc,*vsite_pbc_nalloc;
1185 t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
1186 gmx_reverse_top_t *rt;
1187 gmx_molblock_ind_t *mbi;
1190 if (vsite && vsite->n_intercg_vsite > 0)
1192 vsite_pbc = vsite->vsite_pbc_loc;
1193 vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
1198 vsite_pbc_nalloc = NULL;
1201 /* Clear the counts */
1202 for(ftype=0; ftype<F_NRE; ftype++)
1204 idef->il[ftype].nr = 0;
1208 rt = dd->reverse_top;
1210 if (rt->ril_mt_tot_size == 0)
1212 /* There are no interactions to assign */
1213 return nbonded_local;
1218 for(i=0; i<dd->nat_home; i++)
1220 /* Get the global atom number */
1221 i_gl = dd->gatindex[i];
1222 global_atomnr_to_moltype_ind(mbi,i_gl,&mb,&mt,&mol,&i_mol);
1223 /* Check all interactions assigned to this atom */
1224 index = rt->ril_mt[mt].index;
1225 rtil = rt->ril_mt[mt].il;
1226 /* Check all interactions assigned to this atom */
1228 while (j < index[i_mol+1])
1233 if (interaction_function[ftype].flags & IF_VSITE)
1235 /* The vsite construction goes where the vsite itself is */
1236 add_vsite(dd->ga2la,index,rtil,ftype,nral,
1238 iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
1244 tiatoms[0] = iatoms[0];
1246 for(k=2; k<=nral; k++)
1248 tiatoms[k] = i + iatoms[k] - iatoms[1];
1250 if (ftype == F_POSRES)
1252 add_posres(mol,i_mol,&molb[mb],tiatoms,idef);
1254 else if (ftype == F_FBPOSRES)
1256 add_fbposres(mol,i_mol,&molb[mb],tiatoms,idef);
1258 /* Add this interaction to the local topology */
1259 add_ifunc(nral,tiatoms,&idef->il[ftype]);
1260 /* Sum so we can check in global_stat if we have everything */
1267 return nbonded_local;
1270 static int make_local_exclusions(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1272 gmx_bool bRCheck,real rc,
1273 int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1277 int nizone,n,count,ic,jla0,jla1,jla;
1278 int cg,la0,la1,la,a_gl,mb,mt,mol,a_mol,j,aj_mol;
1283 gmx_molblock_ind_t *mbi;
1286 /* Since for RF and PME we need to loop over the exclusions
1287 * we should store each exclusion only once. This is done
1288 * using the same zone scheme as used for neighbor searching.
1289 * The exclusions involving non-home atoms are stored only
1290 * one way: atom j is in the excl list of i only for j > i,
1291 * where i and j are local atom numbers.
1294 lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1295 if (lexcls->nr+1 > lexcls->nalloc_index)
1297 lexcls->nalloc_index = over_alloc_dd(lexcls->nr)+1;
1298 srenew(lexcls->index,lexcls->nalloc_index);
1301 mbi = dd->reverse_top->mbi;
1307 if (dd->n_intercg_excl)
1309 nizone = zones->nizone;
1317 for(ic=0; ic<nizone; ic++)
1319 jla0 = dd->cgindex[zones->izone[ic].jcg0];
1320 jla1 = dd->cgindex[zones->izone[ic].jcg1];
1321 for(cg=zones->cg_range[ic]; cg<zones->cg_range[ic+1]; cg++)
1323 /* Here we assume the number of exclusions in one charge group
1324 * is never larger than 1000.
1326 if (n+1000 > lexcls->nalloc_a)
1328 lexcls->nalloc_a = over_alloc_large(n+1000);
1329 srenew(lexcls->a,lexcls->nalloc_a);
1331 la0 = dd->cgindex[cg];
1332 la1 = dd->cgindex[cg+1];
1333 if (GET_CGINFO_EXCL_INTER(fr->cginfo[cg]) ||
1334 !GET_CGINFO_EXCL_INTRA(fr->cginfo[cg]))
1336 /* Copy the exclusions from the global top */
1337 for(la=la0; la<la1; la++) {
1338 lexcls->index[la] = n;
1339 a_gl = dd->gatindex[la];
1340 global_atomnr_to_moltype_ind(mbi,a_gl,&mb,&mt,&mol,&a_mol);
1341 excls = &mtop->moltype[mt].excls;
1342 for(j=excls->index[a_mol]; j<excls->index[a_mol+1]; j++)
1344 aj_mol = excls->a[j];
1345 /* This computation of jla is only correct intra-cg */
1346 jla = la + aj_mol - a_mol;
1347 if (jla >= la0 && jla < la1)
1349 /* This is an intra-cg exclusion. We can skip
1350 * the global indexing and distance checking.
1352 /* Intra-cg exclusions are only required
1353 * for the home zone.
1357 lexcls->a[n++] = jla;
1358 /* Check to avoid double counts */
1367 /* This is a inter-cg exclusion */
1368 /* Since exclusions are pair interactions,
1369 * just like non-bonded interactions,
1370 * they can be assigned properly up
1371 * to the DD cutoff (not cutoff_min as
1372 * for the other bonded interactions).
1374 if (ga2la_get(ga2la,a_gl+aj_mol-a_mol,&jla,&cell))
1376 if (ic == 0 && cell == 0)
1378 lexcls->a[n++] = jla;
1379 /* Check to avoid double counts */
1385 else if (jla >= jla0 && jla < jla1 &&
1387 dd_dist2(pbc_null,cg_cm,la2lc,la,jla) < rc2))
1389 /* jla > la, since jla0 > la */
1390 lexcls->a[n++] = jla;
1400 /* There are no inter-cg excls and this cg is self-excluded.
1401 * These exclusions are only required for zone 0,
1402 * since other zones do not see themselves.
1406 for(la=la0; la<la1; la++)
1408 lexcls->index[la] = n;
1409 for(j=la0; j<la1; j++)
1414 count += ((la1 - la0)*(la1 - la0 - 1))/2;
1418 /* We don't need exclusions for this cg */
1419 for(la=la0; la<la1; la++)
1421 lexcls->index[la] = n;
1427 if (dd->n_intercg_excl == 0)
1429 /* There are no exclusions involving non-home charge groups,
1430 * but we need to set the indices for neighborsearching.
1432 la0 = dd->cgindex[zones->izone[0].cg1];
1433 for(la=la0; la<lexcls->nr; la++)
1435 lexcls->index[la] = n;
1438 lexcls->index[lexcls->nr] = n;
1440 if (dd->n_intercg_excl == 0)
1442 /* nr is only used to loop over the exclusions for Ewald and RF,
1443 * so we can set it to the number of home atoms for efficiency.
1445 lexcls->nr = dd->cgindex[zones->izone[0].cg1];
1449 fprintf(debug,"We have %d exclusions, check count %d\n",
1456 void dd_make_local_cgs(gmx_domdec_t *dd,t_block *lcgs)
1458 lcgs->nr = dd->ncg_tot;
1459 lcgs->index = dd->cgindex;
1462 void dd_make_local_top(FILE *fplog,
1463 gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1464 int npbcdim,matrix box,
1465 rvec cellsize_min,ivec npulse,
1466 t_forcerec *fr,gmx_vsite_t *vsite,
1467 gmx_mtop_t *mtop,gmx_localtop_t *ltop)
1469 gmx_bool bUniqueExcl,bRCheckMB,bRCheck2B,bRCheckExcl;
1473 t_pbc pbc,*pbc_null=NULL;
1477 fprintf(debug,"Making local topology\n");
1480 dd_make_local_cgs(dd,<op->cgs);
1484 bRCheckExcl = FALSE;
1486 if (!dd->reverse_top->bMultiCGmols)
1488 /* We don't need checks, assign all interactions with local atoms */
1490 dd->nbonded_local = make_local_bondeds_intracg(dd,mtop->molblock,
1495 /* We need to check to which cell bondeds should be assigned */
1496 rc = dd_cutoff_twobody(dd);
1499 fprintf(debug,"Two-body bonded cut-off distance is %g\n",rc);
1502 /* Should we check cg_cm distances when assigning bonded interactions? */
1503 for(d=0; d<DIM; d++)
1506 /* Only need to check for dimensions where the part of the box
1507 * that is not communicated is smaller than the cut-off.
1509 if (d < npbcdim && dd->nc[d] > 1 &&
1510 (dd->nc[d] - npulse[d])*cellsize_min[d] < 2*rc)
1517 /* Check for interactions between two atoms,
1518 * where we can allow interactions up to the cut-off,
1519 * instead of up to the smallest cell dimension.
1526 "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
1527 d,cellsize_min[d],d,rcheck[d],bRCheck2B);
1530 if (dd->reverse_top->bExclRequired)
1532 bRCheckExcl = bRCheck2B;
1536 /* If we don't have forces on exclusions,
1537 * we don't care about exclusions being assigned mulitple times.
1539 bRCheckExcl = FALSE;
1541 if (bRCheckMB || bRCheck2B)
1546 set_pbc_dd(&pbc,fr->ePBC,dd,TRUE,box);
1555 dd->nbonded_local = make_local_bondeds(dd,zones,mtop->molblock,
1556 bRCheckMB,rcheck,bRCheck2B,rc,
1562 /* The ilist is not sorted yet,
1563 * we can only do this when we have the charge arrays.
1565 ltop->idef.ilsort = ilsortUNKNOWN;
1567 nexcl = make_local_exclusions(dd,zones,mtop,bRCheckExcl,
1568 rc,dd->la2lc,pbc_null,fr->cg_cm,
1571 if (dd->reverse_top->bExclRequired)
1573 dd->nbonded_local += nexcl;
1576 ltop->atomtypes = mtop->atomtypes;
1578 /* For an error message only */
1579 dd->reverse_top->err_top_global = mtop;
1580 dd->reverse_top->err_top_local = ltop;
1583 void dd_sort_local_top(gmx_domdec_t *dd,t_mdatoms *mdatoms,
1584 gmx_localtop_t *ltop)
1586 if (dd->reverse_top->ilsort == ilsortNO_FE)
1588 ltop->idef.ilsort = ilsortNO_FE;
1592 gmx_sort_ilist_fe(<op->idef,mdatoms->chargeA,mdatoms->chargeB);
1596 gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global)
1598 gmx_localtop_t *top;
1603 top->idef.ntypes = top_global->ffparams.ntypes;
1604 top->idef.atnr = top_global->ffparams.atnr;
1605 top->idef.functype = top_global->ffparams.functype;
1606 top->idef.iparams = top_global->ffparams.iparams;
1607 top->idef.fudgeQQ = top_global->ffparams.fudgeQQ;
1608 top->idef.cmap_grid= top_global->ffparams.cmap_grid;
1610 for(i=0; i<F_NRE; i++)
1612 top->idef.il[i].iatoms = NULL;
1613 top->idef.il[i].nalloc = 0;
1615 top->idef.ilsort = ilsortUNKNOWN;
1620 void dd_init_local_state(gmx_domdec_t *dd,
1621 t_state *state_global,t_state *state_local)
1623 int buf[NITEM_DD_INIT_LOCAL_STATE];
1627 buf[0] = state_global->flags;
1628 buf[1] = state_global->ngtc;
1629 buf[2] = state_global->nnhpres;
1630 buf[3] = state_global->nhchainlength;
1631 buf[4] = state_global->dfhist.nlambda;
1633 dd_bcast(dd,NITEM_DD_INIT_LOCAL_STATE*sizeof(int),buf);
1635 init_state(state_local,0,buf[1],buf[2],buf[3],buf[4]);
1636 state_local->flags = buf[0];
1638 /* With Langevin Dynamics we need to make proper storage space
1639 * in the global and local state for the random numbers.
1641 if (state_local->flags & (1<<estLD_RNG))
1643 if (DDMASTER(dd) && state_global->nrngi > 1)
1645 state_global->nrng = dd->nnodes*gmx_rng_n();
1646 srenew(state_global->ld_rng,state_global->nrng);
1648 state_local->nrng = gmx_rng_n();
1649 snew(state_local->ld_rng,state_local->nrng);
1651 if (state_local->flags & (1<<estLD_RNGI))
1653 if (DDMASTER(dd) && state_global->nrngi > 1)
1655 state_global->nrngi = dd->nnodes;
1656 srenew(state_global->ld_rngi,state_global->nrngi);
1658 snew(state_local->ld_rngi,1);
1662 static void check_link(t_blocka *link,int cg_gl,int cg_gl_j)
1668 for(k=link->index[cg_gl]; k<link->index[cg_gl+1]; k++)
1670 if (link->a[k] == cg_gl_j)
1677 /* Add this charge group link */
1678 if (link->index[cg_gl+1]+1 > link->nalloc_a)
1680 link->nalloc_a = over_alloc_large(link->index[cg_gl+1]+1);
1681 srenew(link->a,link->nalloc_a);
1683 link->a[link->index[cg_gl+1]] = cg_gl_j;
1684 link->index[cg_gl+1]++;
1688 static int *make_at2cg(t_block *cgs)
1692 snew(at2cg,cgs->index[cgs->nr]);
1693 for(cg=0; cg<cgs->nr; cg++)
1695 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1704 t_blocka *make_charge_group_links(gmx_mtop_t *mtop,gmx_domdec_t *dd,
1705 cginfo_mb_t *cginfo_mb)
1707 gmx_reverse_top_t *rt;
1708 int mb,cg_offset,cg,cg_gl,a,aj,i,j,ftype,nral,nlink_mol,mol,ncgi;
1709 gmx_molblock_t *molb;
1710 gmx_moltype_t *molt;
1714 gmx_reverse_ilist_t ril;
1716 cginfo_mb_t *cgi_mb;
1718 /* For each charge group make a list of other charge groups
1719 * in the system that a linked to it via bonded interactions
1720 * which are also stored in reverse_top.
1723 rt = dd->reverse_top;
1726 snew(link->index,ncg_mtop(mtop)+1);
1733 for(mb=0; mb<mtop->nmolblock; mb++)
1735 molb = &mtop->molblock[mb];
1736 if (molb->nmol == 0)
1740 molt = &mtop->moltype[molb->type];
1742 excls = &molt->excls;
1743 a2c = make_at2cg(cgs);
1744 /* Make a reverse ilist in which the interactions are linked
1745 * to all atoms, not only the first atom as in gmx_reverse_top.
1746 * The constraints are discarded here.
1748 make_reverse_ilist(molt,NULL,FALSE,FALSE,TRUE,&ril);
1750 cgi_mb = &cginfo_mb[mb];
1752 for(cg=0; cg<cgs->nr; cg++)
1754 cg_gl = cg_offset + cg;
1755 link->index[cg_gl+1] = link->index[cg_gl];
1756 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1759 while (i < ril.index[a+1])
1761 ftype = ril.il[i++];
1763 /* Skip the ifunc index */
1765 for(j=0; j<nral; j++)
1770 check_link(link,cg_gl,cg_offset+a2c[aj]);
1773 i += nral_rt(ftype);
1775 if (rt->bExclRequired)
1777 /* Exclusions always go both ways */
1778 for(j=excls->index[a]; j<excls->index[a+1]; j++)
1783 check_link(link,cg_gl,cg_offset+a2c[aj]);
1788 if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
1790 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
1794 nlink_mol = link->index[cg_offset+cgs->nr] - link->index[cg_offset];
1796 cg_offset += cgs->nr;
1798 destroy_reverse_ilist(&ril);
1803 fprintf(debug,"molecule type '%s' %d cgs has %d cg links through bonded interac.\n",*molt->name,cgs->nr,nlink_mol);
1808 /* Copy the data for the rest of the molecules in this block */
1809 link->nalloc_a += (molb->nmol - 1)*nlink_mol;
1810 srenew(link->a,link->nalloc_a);
1811 for(mol=1; mol<molb->nmol; mol++)
1813 for(cg=0; cg<cgs->nr; cg++)
1815 cg_gl = cg_offset + cg;
1816 link->index[cg_gl+1] =
1817 link->index[cg_gl+1-cgs->nr] + nlink_mol;
1818 for(j=link->index[cg_gl]; j<link->index[cg_gl+1]; j++)
1820 link->a[j] = link->a[j-nlink_mol] + cgs->nr;
1822 if (link->index[cg_gl+1] - link->index[cg_gl] > 0 &&
1823 cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
1825 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
1829 cg_offset += cgs->nr;
1836 fprintf(debug,"Of the %d charge groups %d are linked via bonded interactions\n",ncg_mtop(mtop),ncgi);
1842 static void bonded_cg_distance_mol(gmx_moltype_t *molt,int *at2cg,
1843 gmx_bool bBCheck,gmx_bool bExcl,rvec *cg_cm,
1844 real *r_2b,int *ft2b,int *a2_1,int *a2_2,
1845 real *r_mb,int *ftmb,int *am_1,int *am_2)
1847 int ftype,nral,i,j,ai,aj,cgi,cgj;
1850 real r2_2b,r2_mb,rij2;
1854 for(ftype=0; ftype<F_NRE; ftype++)
1856 if (dd_check_ftype(ftype,bBCheck,FALSE))
1858 il = &molt->ilist[ftype];
1862 for(i=0; i<il->nr; i+=1+nral)
1864 for(ai=0; ai<nral; ai++)
1866 cgi = at2cg[il->iatoms[i+1+ai]];
1867 for(aj=0; aj<nral; aj++) {
1868 cgj = at2cg[il->iatoms[i+1+aj]];
1871 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
1872 if (nral == 2 && rij2 > r2_2b)
1876 *a2_1 = il->iatoms[i+1+ai];
1877 *a2_2 = il->iatoms[i+1+aj];
1879 if (nral > 2 && rij2 > r2_mb)
1883 *am_1 = il->iatoms[i+1+ai];
1884 *am_2 = il->iatoms[i+1+aj];
1895 excls = &molt->excls;
1896 for(ai=0; ai<excls->nr; ai++)
1899 for(j=excls->index[ai]; j<excls->index[ai+1]; j++) {
1900 cgj = at2cg[excls->a[j]];
1903 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
1913 *r_2b = sqrt(r2_2b);
1914 *r_mb = sqrt(r2_mb);
1917 static void get_cgcm_mol(gmx_moltype_t *molt,gmx_ffparams_t *ffparams,
1918 int ePBC,t_graph *graph,matrix box,
1920 rvec *x,rvec *xs,rvec *cg_cm)
1924 if (ePBC != epbcNONE)
1926 mk_mshift(NULL,graph,ePBC,box,x);
1928 shift_x(graph,box,x,xs);
1929 /* By doing an extra mk_mshift the molecules that are broken
1930 * because they were e.g. imported from another software
1931 * will be made whole again. Such are the healing powers
1934 mk_mshift(NULL,graph,ePBC,box,xs);
1938 /* We copy the coordinates so the original coordinates remain
1939 * unchanged, just to be 100% sure that we do not affect
1940 * binary reproducibility of simulations.
1942 n = molt->cgs.index[molt->cgs.nr];
1945 copy_rvec(x[i],xs[i]);
1951 construct_vsites(NULL,vsite,xs,NULL,0.0,NULL,
1952 ffparams->iparams,molt->ilist,
1953 epbcNONE,TRUE,NULL,NULL,NULL);
1956 calc_cgcm(NULL,0,molt->cgs.nr,&molt->cgs,xs,cg_cm);
1959 static int have_vsite_molt(gmx_moltype_t *molt)
1965 for(i=0; i<F_NRE; i++)
1967 if ((interaction_function[i].flags & IF_VSITE) &&
1968 molt->ilist[i].nr > 0) {
1976 void dd_bonded_cg_distance(FILE *fplog,
1977 gmx_domdec_t *dd,gmx_mtop_t *mtop,
1978 t_inputrec *ir,rvec *x,matrix box,
1980 real *r_2b,real *r_mb)
1982 gmx_bool bExclRequired;
1983 int mb,cg_offset,at_offset,*at2cg,mol;
1986 gmx_molblock_t *molb;
1987 gmx_moltype_t *molt;
1989 real rmol_2b,rmol_mb;
1990 int ft2b=-1,a_2b_1=-1,a_2b_2=-1,ftmb=-1,a_mb_1=-1,a_mb_2=-1;
1991 int ftm2b=-1,amol_2b_1=-1,amol_2b_2=-1,ftmmb=-1,amol_mb_1=-1,amol_mb_2=-1;
1993 bExclRequired = IR_EXCL_FORCES(*ir);
1995 /* For gmx_vsite_t everything 0 should work (without pbc) */
2002 for(mb=0; mb<mtop->nmolblock; mb++)
2004 molb = &mtop->molblock[mb];
2005 molt = &mtop->moltype[molb->type];
2006 if (molt->cgs.nr == 1 || molb->nmol == 0)
2008 cg_offset += molb->nmol*molt->cgs.nr;
2009 at_offset += molb->nmol*molt->atoms.nr;
2013 if (ir->ePBC != epbcNONE)
2015 mk_graph_ilist(NULL,molt->ilist,0,molt->atoms.nr,FALSE,FALSE,
2019 at2cg = make_at2cg(&molt->cgs);
2020 snew(xs,molt->atoms.nr);
2021 snew(cg_cm,molt->cgs.nr);
2022 for(mol=0; mol<molb->nmol; mol++)
2024 get_cgcm_mol(molt,&mtop->ffparams,ir->ePBC,&graph,box,
2025 have_vsite_molt(molt) ? vsite : NULL,
2026 x+at_offset,xs,cg_cm);
2028 bonded_cg_distance_mol(molt,at2cg,bBCheck,bExclRequired,cg_cm,
2029 &rmol_2b,&ftm2b,&amol_2b_1,&amol_2b_2,
2030 &rmol_mb,&ftmmb,&amol_mb_1,&amol_mb_2);
2031 if (rmol_2b > *r_2b)
2035 a_2b_1 = at_offset + amol_2b_1;
2036 a_2b_2 = at_offset + amol_2b_2;
2038 if (rmol_mb > *r_mb)
2042 a_mb_1 = at_offset + amol_mb_1;
2043 a_mb_2 = at_offset + amol_mb_2;
2046 cg_offset += molt->cgs.nr;
2047 at_offset += molt->atoms.nr;
2052 if (ir->ePBC != epbcNONE)
2061 if (fplog && (ft2b >= 0 || ftmb >= 0))
2064 "Initial maximum inter charge-group distances:\n");
2068 " two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2069 *r_2b,interaction_function[ft2b].longname,
2075 " multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2076 *r_mb,interaction_function[ftmb].longname,