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 #include "gmx_omp_nthreads.h"
42 /* for dd_init_local_state */
43 #define NITEM_DD_INIT_LOCAL_STATE 5
46 int *index; /* Index for each atom into il */
47 int *il; /* ftype|type|a0|...|an|ftype|... */
48 } gmx_reverse_ilist_t;
57 typedef struct gmx_reverse_top {
58 gmx_bool bExclRequired; /* Do we require all exclusions to be assigned? */
59 gmx_bool bConstr; /* Are there constraints in this revserse top? */
60 gmx_bool bSettle; /* Are there settles in this revserse top? */
61 gmx_bool bBCheck; /* All bonded interactions have to be assigned? */
62 gmx_bool bMultiCGmols; /* Are the multi charge-group molecules? */
63 gmx_reverse_ilist_t *ril_mt; /* Reverse ilist for all moltypes */
65 int ilsort; /* The sorting state of bondeds for free energy */
66 gmx_molblock_ind_t *mbi;
69 /* Work data structures for multi-threading */
73 int **vsite_pbc_nalloc;
75 t_blocka *excl_thread;
76 int *excl_count_thread;
78 /* Pointers only used for an error message */
79 gmx_mtop_t *err_top_global;
80 gmx_localtop_t *err_top_local;
83 static int nral_rt(int ftype)
85 /* Returns the number of atom entries for il in gmx_reverse_top_t */
89 if (interaction_function[ftype].flags & IF_VSITE)
91 /* With vsites the reverse topology contains
92 * two extra entries for PBC.
100 /* This function tells which interactions need to be assigned exactly once */
101 static gmx_bool dd_check_ftype(int ftype,gmx_bool bBCheck,
102 gmx_bool bConstr,gmx_bool bSettle)
104 return (((interaction_function[ftype].flags & IF_BOND) &&
105 !(interaction_function[ftype].flags & IF_VSITE) &&
106 (bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO))) ||
107 (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) ||
108 (bSettle && ftype == F_SETTLE));
111 static void print_error_header(FILE *fplog,char *moltypename,int nprint)
113 fprintf(fplog, "\nMolecule type '%s'\n",moltypename);
114 fprintf(stderr,"\nMolecule type '%s'\n",moltypename);
116 "the first %d missing interactions, except for exclusions:\n",
119 "the first %d missing interactions, except for exclusions:\n",
123 static void print_missing_interactions_mb(FILE *fplog,t_commrec *cr,
124 gmx_reverse_top_t *rt,
126 gmx_reverse_ilist_t *ril,
127 int a_start,int a_end,
128 int nat_mol,int nmol,
131 int nril_mol,*assigned,*gatindex;
132 int ftype,ftype_j,nral,i,j_mol,j,k,a0,a0_mol,mol,a,a_gl;
138 nril_mol = ril->index[nat_mol];
139 snew(assigned,nmol*nril_mol);
141 gatindex = cr->dd->gatindex;
142 for(ftype=0; ftype<F_NRE; ftype++)
144 if (dd_check_ftype(ftype,rt->bBCheck,rt->bConstr,rt->bSettle))
147 il = &idef->il[ftype];
149 for(i=0; i<il->nr; i+=1+nral)
151 a0 = gatindex[ia[1]];
152 /* Check if this interaction is in
153 * the currently checked molblock.
155 if (a0 >= a_start && a0 < a_end)
157 mol = (a0 - a_start)/nat_mol;
158 a0_mol = (a0 - a_start) - mol*nat_mol;
159 j_mol = ril->index[a0_mol];
161 while (j_mol < ril->index[a0_mol+1] && !bFound)
163 j = mol*nril_mol + j_mol;
164 ftype_j = ril->il[j_mol];
165 /* Here we need to check if this interaction has
166 * not already been assigned, since we could have
167 * multiply defined interactions.
169 if (ftype == ftype_j && ia[0] == ril->il[j_mol+1] &&
172 /* Check the atoms */
174 for(a=0; a<nral; a++)
176 if (gatindex[ia[1+a]] !=
177 a_start + mol*nat_mol + ril->il[j_mol+2+a])
187 j_mol += 2 + nral_rt(ftype_j);
191 gmx_incons("Some interactions seem to be assigned multiple times");
199 gmx_sumi(nmol*nril_mol,assigned,cr);
203 for(mol=0; mol<nmol; mol++)
206 while (j_mol < nril_mol)
208 ftype = ril->il[j_mol];
210 j = mol*nril_mol + j_mol;
211 if (assigned[j] == 0 &&
212 !(interaction_function[ftype].flags & IF_VSITE))
214 if (DDMASTER(cr->dd))
218 print_error_header(fplog,moltypename,nprint);
220 fprintf(fplog, "%20s atoms",
221 interaction_function[ftype].longname);
222 fprintf(stderr,"%20s atoms",
223 interaction_function[ftype].longname);
224 for(a=0; a<nral; a++) {
225 fprintf(fplog, "%5d",ril->il[j_mol+2+a]+1);
226 fprintf(stderr,"%5d",ril->il[j_mol+2+a]+1);
234 fprintf(fplog, " global");
235 fprintf(stderr," global");
236 for(a=0; a<nral; a++)
238 fprintf(fplog, "%6d",
239 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
240 fprintf(stderr,"%6d",
241 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
243 fprintf(fplog, "\n");
244 fprintf(stderr,"\n");
252 j_mol += 2 + nral_rt(ftype);
259 static void print_missing_interactions_atoms(FILE *fplog,t_commrec *cr,
260 gmx_mtop_t *mtop,t_idef *idef)
262 int mb,a_start,a_end;
263 gmx_molblock_t *molb;
264 gmx_reverse_top_t *rt;
266 rt = cr->dd->reverse_top;
268 /* Print the atoms in the missing interactions per molblock */
270 for(mb=0; mb<mtop->nmolblock; mb++)
272 molb = &mtop->molblock[mb];
274 a_end = a_start + molb->nmol*molb->natoms_mol;
276 print_missing_interactions_mb(fplog,cr,rt,
277 *(mtop->moltype[molb->type].name),
278 &rt->ril_mt[molb->type],
279 a_start,a_end,molb->natoms_mol,
285 void dd_print_missing_interactions(FILE *fplog,t_commrec *cr,int local_count, gmx_mtop_t *top_global, t_state *state_local)
287 int ndiff_tot,cl[F_NRE],n,ndiff,rest_global,rest_local;
291 gmx_mtop_t *err_top_global;
292 gmx_localtop_t *err_top_local;
296 err_top_global = dd->reverse_top->err_top_global;
297 err_top_local = dd->reverse_top->err_top_local;
301 fprintf(fplog,"\nNot all bonded interactions have been properly assigned to the domain decomposition cells\n");
305 ndiff_tot = local_count - dd->nbonded_global;
307 for(ftype=0; ftype<F_NRE; ftype++)
310 cl[ftype] = err_top_local->idef.il[ftype].nr/(1+nral);
313 gmx_sumi(F_NRE,cl,cr);
317 fprintf(fplog,"\nA list of missing interactions:\n");
318 fprintf(stderr,"\nA list of missing interactions:\n");
319 rest_global = dd->nbonded_global;
320 rest_local = local_count;
321 for(ftype=0; ftype<F_NRE; ftype++)
323 /* In the reverse and local top all constraints are merged
324 * into F_CONSTR. So in the if statement we skip F_CONSTRNC
325 * and add these constraints when doing F_CONSTR.
327 if (((interaction_function[ftype].flags & IF_BOND) &&
328 (dd->reverse_top->bBCheck
329 || !(interaction_function[ftype].flags & IF_LIMZERO)))
330 || (dd->reverse_top->bConstr && ftype == F_CONSTR)
331 || (dd->reverse_top->bSettle && ftype == F_SETTLE))
334 n = gmx_mtop_ftype_count(err_top_global,ftype);
335 if (ftype == F_CONSTR)
337 n += gmx_mtop_ftype_count(err_top_global,F_CONSTRNC);
339 ndiff = cl[ftype] - n;
342 sprintf(buf,"%20s of %6d missing %6d",
343 interaction_function[ftype].longname,n,-ndiff);
344 fprintf(fplog,"%s\n",buf);
345 fprintf(stderr,"%s\n",buf);
348 rest_local -= cl[ftype];
352 ndiff = rest_local - rest_global;
355 sprintf(buf,"%20s of %6d missing %6d","exclusions",
357 fprintf(fplog,"%s\n",buf);
358 fprintf(stderr,"%s\n",buf);
362 print_missing_interactions_atoms(fplog,cr,err_top_global,
363 &err_top_local->idef);
364 write_dd_pdb("dd_dump_err",0,"dump",top_global,cr,
365 -1,state_local->x,state_local->box);
370 gmx_incons("One or more interactions were multiple assigned in the domain decompostion");
374 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));
379 static void global_atomnr_to_moltype_ind(gmx_reverse_top_t *rt,int i_gl,
380 int *mb,int *mt,int *mol,int *i_mol)
385 gmx_molblock_ind_t *mbi = rt->mbi;
387 int end = rt->nmolblock; /* exclusive */
390 /* binary search for molblock_ind */
393 if (i_gl >= mbi[mid].a_end)
397 else if (i_gl < mbi[mid].a_start)
411 *mol = (i_gl - mbi->a_start) / mbi->natoms_mol;
412 *i_mol = (i_gl - mbi->a_start) - (*mol)*mbi->natoms_mol;
415 static int count_excls(t_block *cgs,t_blocka *excls,int *n_intercg_excl)
417 int n,n_inter,cg,at0,at1,at,excl,atj;
421 for(cg=0; cg<cgs->nr; cg++)
423 at0 = cgs->index[cg];
424 at1 = cgs->index[cg+1];
425 for(at=at0; at<at1; at++)
427 for(excl=excls->index[at]; excl<excls->index[at+1]; excl++)
429 atj = excls->a[excl];
433 if (atj < at0 || atj >= at1)
445 static int low_make_reverse_ilist(t_ilist *il_mt,t_atom *atom,
448 gmx_bool bConstr,gmx_bool bSettle,
450 int *r_index,int *r_il,
451 gmx_bool bLinkToAllAtoms,
454 int ftype,nral,i,j,nlink,link;
462 for(ftype=0; ftype<F_NRE; ftype++)
464 if ((interaction_function[ftype].flags & (IF_BOND | IF_VSITE)) ||
465 (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) ||
466 (bSettle && ftype == F_SETTLE))
468 bVSite = (interaction_function[ftype].flags & IF_VSITE);
472 for(i=0; i<il->nr; i+=1+nral)
479 /* We don't need the virtual sites for the cg-links */
489 /* Couple to the first atom in the interaction */
492 for(link=0; link<nlink; link++)
497 r_il[r_index[a]+count[a]] =
498 (ftype == F_CONSTRNC ? F_CONSTR : ftype);
499 r_il[r_index[a]+count[a]+1] = ia[0];
500 for(j=1; j<1+nral; j++)
502 /* Store the molecular atom number */
503 r_il[r_index[a]+count[a]+1+j] = ia[j];
506 if (interaction_function[ftype].flags & IF_VSITE)
510 /* Add an entry to iatoms for storing
511 * which of the constructing atoms are
514 r_il[r_index[a]+count[a]+2+nral] = 0;
515 for(j=2; j<1+nral; j++)
517 if (atom[ia[j]].ptype == eptVSite)
519 r_il[r_index[a]+count[a]+2+nral] |= (2<<j);
522 /* Store vsite pbc atom in a second extra entry */
523 r_il[r_index[a]+count[a]+2+nral+1] =
524 (vsite_pbc ? vsite_pbc[ftype-F_VSITE2][i/(1+nral)] : -2);
529 /* We do not count vsites since they are always
530 * uniquely assigned and can be assigned
531 * to multiple nodes with recursive vsites.
534 !(interaction_function[ftype].flags & IF_LIMZERO))
539 count[a] += 2 + nral_rt(ftype);
548 static int make_reverse_ilist(gmx_moltype_t *molt,
550 gmx_bool bConstr,gmx_bool bSettle,
552 gmx_bool bLinkToAllAtoms,
553 gmx_reverse_ilist_t *ril_mt)
555 int nat_mt,*count,i,nint_mt;
557 /* Count the interactions */
558 nat_mt = molt->atoms.nr;
560 low_make_reverse_ilist(molt->ilist,molt->atoms.atom,vsite_pbc,
562 bConstr,bSettle,bBCheck,NULL,NULL,
563 bLinkToAllAtoms,FALSE);
565 snew(ril_mt->index,nat_mt+1);
566 ril_mt->index[0] = 0;
567 for(i=0; i<nat_mt; i++)
569 ril_mt->index[i+1] = ril_mt->index[i] + count[i];
572 snew(ril_mt->il,ril_mt->index[nat_mt]);
574 /* Store the interactions */
576 low_make_reverse_ilist(molt->ilist,molt->atoms.atom,vsite_pbc,
578 bConstr,bSettle,bBCheck,
579 ril_mt->index,ril_mt->il,
580 bLinkToAllAtoms,TRUE);
587 static void destroy_reverse_ilist(gmx_reverse_ilist_t *ril)
593 static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop,gmx_bool bFE,
594 int ***vsite_pbc_molt,
595 gmx_bool bConstr,gmx_bool bSettle,
596 gmx_bool bBCheck,int *nint)
599 gmx_reverse_top_t *rt;
606 /* Should we include constraints (for SHAKE) in rt? */
607 rt->bConstr = bConstr;
608 rt->bSettle = bSettle;
609 rt->bBCheck = bBCheck;
611 rt->bMultiCGmols = FALSE;
612 snew(nint_mt,mtop->nmoltype);
613 snew(rt->ril_mt,mtop->nmoltype);
614 rt->ril_mt_tot_size = 0;
615 for(mt=0; mt<mtop->nmoltype; mt++)
617 molt = &mtop->moltype[mt];
618 if (molt->cgs.nr > 1)
620 rt->bMultiCGmols = TRUE;
623 /* Make the atom to interaction list for this molecule type */
625 make_reverse_ilist(molt,vsite_pbc_molt ? vsite_pbc_molt[mt] : NULL,
626 rt->bConstr,rt->bSettle,rt->bBCheck,FALSE,
629 rt->ril_mt_tot_size += rt->ril_mt[mt].index[molt->atoms.nr];
633 fprintf(debug,"The total size of the atom to interaction index is %d integers\n",rt->ril_mt_tot_size);
637 for(mb=0; mb<mtop->nmolblock; mb++)
639 *nint += mtop->molblock[mb].nmol*nint_mt[mtop->molblock[mb].type];
643 if (bFE && gmx_mtop_bondeds_free_energy(mtop))
645 rt->ilsort = ilsortFE_UNSORTED;
648 rt->ilsort = ilsortNO_FE;
651 /* Make a molblock index for fast searching */
652 snew(rt->mbi,mtop->nmolblock);
653 rt->nmolblock = mtop->nmolblock;
655 for(mb=0; mb<mtop->nmolblock; mb++)
657 rt->mbi[mb].a_start = i;
658 i += mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
659 rt->mbi[mb].a_end = i;
660 rt->mbi[mb].natoms_mol = mtop->molblock[mb].natoms_mol;
661 rt->mbi[mb].type = mtop->molblock[mb].type;
664 rt->nthread = gmx_omp_nthreads_get(emntDomdec);
665 snew(rt->idef_thread,rt->nthread);
666 if (vsite_pbc_molt != NULL)
668 snew(rt->vsite_pbc,rt->nthread);
669 snew(rt->vsite_pbc_nalloc,rt->nthread);
670 for(thread=0; thread<rt->nthread; thread++)
672 snew(rt->vsite_pbc[thread],F_VSITEN-F_VSITE2+1);
673 snew(rt->vsite_pbc_nalloc[thread],F_VSITEN-F_VSITE2+1);
676 snew(rt->nbonded_thread,rt->nthread);
677 snew(rt->excl_thread,rt->nthread);
678 snew(rt->excl_count_thread,rt->nthread);
683 void dd_make_reverse_top(FILE *fplog,
684 gmx_domdec_t *dd,gmx_mtop_t *mtop,
685 gmx_vsite_t *vsite,gmx_constr_t constr,
686 t_inputrec *ir,gmx_bool bBCheck)
688 int mb,n_recursive_vsite,nexcl,nexcl_icg,a;
689 gmx_molblock_t *molb;
694 fprintf(fplog,"\nLinking all bonded interactions to atoms\n");
697 /* If normal and/or settle constraints act only within charge groups,
698 * we can store them in the reverse top and simply assign them to domains.
699 * Otherwise we need to assign them to multiple domains and set up
700 * the parallel version constraint algoirthm(s).
703 dd->reverse_top = make_reverse_top(mtop,ir->efep!=efepNO,
704 vsite ? vsite->vsite_pbc_molt : NULL,
705 !dd->bInterCGcons,!dd->bInterCGsettles,
706 bBCheck,&dd->nbonded_global);
708 if (dd->reverse_top->ril_mt_tot_size >= 200000 &&
710 mtop->nmolblock == 1 && mtop->molblock[0].nmol == 1)
712 /* mtop comes from a pre Gromacs 4 tpr file */
713 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";
716 fprintf(fplog,"\n%s\n\n",note);
720 fprintf(stderr,"\n%s\n\n",note);
724 dd->reverse_top->bExclRequired = IR_EXCL_FORCES(*ir);
727 dd->n_intercg_excl = 0;
728 for(mb=0; mb<mtop->nmolblock; mb++)
730 molb = &mtop->molblock[mb];
731 molt = &mtop->moltype[molb->type];
732 nexcl += molb->nmol*count_excls(&molt->cgs,&molt->excls,&nexcl_icg);
733 dd->n_intercg_excl += molb->nmol*nexcl_icg;
735 if (dd->reverse_top->bExclRequired)
737 dd->nbonded_global += nexcl;
738 if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0 && fplog)
740 fprintf(fplog,"There are %d inter charge-group exclusions,\n"
741 "will use an extra communication step for exclusion forces for %s\n",
742 dd->n_intercg_excl,eel_names[ir->coulombtype]);
746 if (vsite && vsite->n_intercg_vsite > 0)
750 fprintf(fplog,"There are %d inter charge-group virtual sites,\n"
751 "will an extra communication step for selected coordinates and forces\n",
752 vsite->n_intercg_vsite);
754 init_domdec_vsites(dd,vsite->n_intercg_vsite);
757 if (dd->bInterCGcons || dd->bInterCGsettles)
759 init_domdec_constraints(dd,mtop,constr);
767 static inline void add_ifunc(int nral,t_iatom *tiatoms,t_ilist *il)
772 if (il->nr+1+nral > il->nalloc)
774 il->nalloc = over_alloc_large(il->nr+1+nral);
775 srenew(il->iatoms,il->nalloc);
777 liatoms = il->iatoms + il->nr;
778 for(k=0; k<=nral; k++)
780 liatoms[k] = tiatoms[k];
785 static void add_posres(int mol,int a_mol,const gmx_molblock_t *molb,
786 t_iatom *iatoms,const t_iparams *ip_in,
792 /* This position restraint has not been added yet,
793 * so it's index is the current number of position restraints.
795 n = idef->il[F_POSRES].nr/2;
796 if (n+1 > idef->iparams_posres_nalloc)
798 idef->iparams_posres_nalloc = over_alloc_dd(n+1);
799 srenew(idef->iparams_posres,idef->iparams_posres_nalloc);
801 ip = &idef->iparams_posres[n];
802 /* Copy the force constants */
803 *ip = ip_in[iatoms[0]];
805 /* Get the position restraint coordinates from the molblock */
806 a_molb = mol*molb->natoms_mol + a_mol;
807 if (a_molb >= molb->nposres_xA)
809 gmx_incons("Not enough position restraint coordinates");
811 ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
812 ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
813 ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
814 if (molb->nposres_xB > 0)
816 ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
817 ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
818 ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
822 ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
823 ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
824 ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
826 /* Set the parameter index for idef->iparams_posre */
830 static void add_fbposres(int mol,int a_mol,const gmx_molblock_t *molb,
831 t_iatom *iatoms,const t_iparams *ip_in,
837 /* This flat-bottom position restraint has not been added yet,
838 * so it's index is the current number of position restraints.
840 n = idef->il[F_FBPOSRES].nr/2;
841 if (n+1 > idef->iparams_fbposres_nalloc)
843 idef->iparams_fbposres_nalloc = over_alloc_dd(n+1);
844 srenew(idef->iparams_fbposres,idef->iparams_fbposres_nalloc);
846 ip = &idef->iparams_fbposres[n];
847 /* Copy the force constants */
848 *ip = ip_in[iatoms[0]];
850 /* Get the position restriant coordinats from the molblock */
851 a_molb = mol*molb->natoms_mol + a_mol;
852 if (a_molb >= molb->nposres_xA)
854 gmx_incons("Not enough position restraint coordinates");
856 /* Take reference positions from A position of normal posres */
857 ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
858 ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
859 ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
861 /* Note: no B-type for flat-bottom posres */
863 /* Set the parameter index for idef->iparams_posre */
867 static void add_vsite(gmx_ga2la_t ga2la,int *index,int *rtil,
869 gmx_bool bHomeA,int a,int a_gl,int a_mol,
871 t_idef *idef,int **vsite_pbc,int *vsite_pbc_nalloc)
873 int k,ak_gl,vsi,pbc_a_mol;
874 t_iatom tiatoms[1+MAXATOMLIST],*iatoms_r;
875 int j,ftype_r,nral_r;
878 tiatoms[0] = iatoms[0];
882 /* We know the local index of the first atom */
887 /* Convert later in make_local_vsites */
888 tiatoms[1] = -a_gl - 1;
891 for(k=2; k<1+nral; k++)
893 ak_gl = a_gl + iatoms[k] - a_mol;
894 if (!ga2la_get_home(ga2la,ak_gl,&tiatoms[k]))
896 /* Copy the global index, convert later in make_local_vsites */
897 tiatoms[k] = -(ak_gl + 1);
901 /* Add this interaction to the local topology */
902 add_ifunc(nral,tiatoms,&idef->il[ftype]);
905 vsi = idef->il[ftype].nr/(1+nral) - 1;
906 if (vsi >= vsite_pbc_nalloc[ftype-F_VSITE2])
908 vsite_pbc_nalloc[ftype-F_VSITE2] = over_alloc_large(vsi+1);
909 srenew(vsite_pbc[ftype-F_VSITE2],vsite_pbc_nalloc[ftype-F_VSITE2]);
913 pbc_a_mol = iatoms[1+nral+1];
916 /* The pbc flag is one of the following two options:
917 * -2: vsite and all constructing atoms are within the same cg, no pbc
918 * -1: vsite and its first constructing atom are in the same cg, do pbc
920 vsite_pbc[ftype-F_VSITE2][vsi] = pbc_a_mol;
924 /* Set the pbc atom for this vsite so we can make its pbc
925 * identical to the rest of the atoms in its charge group.
926 * Since the order of the atoms does not change within a charge
927 * group, we do not need the global to local atom index.
929 vsite_pbc[ftype-F_VSITE2][vsi] = a + pbc_a_mol - iatoms[1];
934 /* This vsite is non-home (required for recursion),
935 * and therefore there is no charge group to match pbc with.
936 * But we always turn on full_pbc to assure that higher order
937 * recursion works correctly.
939 vsite_pbc[ftype-F_VSITE2][vsi] = -1;
945 /* Check for recursion */
946 for(k=2; k<1+nral; k++)
948 if ((iatoms[1+nral] & (2<<k)) && (tiatoms[k] < 0))
950 /* This construction atoms is a vsite and not a home atom */
953 fprintf(debug,"Constructing atom %d of vsite atom %d is a vsite and non-home\n",iatoms[k]+1,a_mol+1);
955 /* Find the vsite construction */
957 /* Check all interactions assigned to this atom */
958 j = index[iatoms[k]];
959 while (j < index[iatoms[k]+1])
962 nral_r = NRAL(ftype_r);
963 if (interaction_function[ftype_r].flags & IF_VSITE)
965 /* Add this vsite (recursion) */
966 add_vsite(ga2la,index,rtil,ftype_r,nral_r,
967 FALSE,-1,a_gl+iatoms[k]-iatoms[1],iatoms[k],
968 rtil+j,idef,vsite_pbc,vsite_pbc_nalloc);
981 static void make_la2lc(gmx_domdec_t *dd)
983 int *cgindex,*la2lc,cg,a;
985 cgindex = dd->cgindex;
987 if (dd->nat_tot > dd->la2lc_nalloc)
989 dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
990 snew(dd->la2lc,dd->la2lc_nalloc);
994 /* Make the local atom to local cg index */
995 for(cg=0; cg<dd->ncg_tot; cg++)
997 for(a=cgindex[cg]; a<cgindex[cg+1]; a++)
1004 static real dd_dist2(t_pbc *pbc_null,rvec *cg_cm,const int *la2lc,int i,int j)
1010 pbc_dx_aiuc(pbc_null,cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
1014 rvec_sub(cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
1020 /* Append the nsrc t_blocka block structures in src to *dest */
1021 static void combine_blocka(t_blocka *dest,const t_blocka *src,int nsrc)
1025 ni = src[nsrc-1].nr;
1027 for(s=0; s<nsrc; s++)
1031 if (ni + 1 > dest->nalloc_index)
1033 dest->nalloc_index = over_alloc_large(ni+1);
1034 srenew(dest->index,dest->nalloc_index);
1036 if (dest->nra + na > dest->nalloc_a)
1038 dest->nalloc_a = over_alloc_large(dest->nra+na);
1039 srenew(dest->a,dest->nalloc_a);
1041 for(s=0; s<nsrc; s++)
1043 for(i=dest->nr+1; i<src[s].nr+1; i++)
1045 dest->index[i] = dest->nra + src[s].index[i];
1047 for(i=0; i<src[s].nra; i++)
1049 dest->a[dest->nra+i] = src[s].a[i];
1051 dest->nr = src[s].nr;
1052 dest->nra += src[s].nra;
1056 /* Append the nsrc t_idef structures in src to *dest,
1057 * virtual sites need special attention, as pbc info differs per vsite.
1059 static void combine_idef(t_idef *dest,const t_idef *src,int nsrc,
1060 gmx_vsite_t *vsite,int ***vsite_pbc_t)
1068 for(ftype=0; ftype<F_NRE; ftype++)
1071 for(s=0; s<nsrc; s++)
1073 n += src[s].il[ftype].nr;
1077 ild = &dest->il[ftype];
1079 if (ild->nr + n > ild->nalloc)
1081 ild->nalloc = over_alloc_large(ild->nr+n);
1082 srenew(ild->iatoms,ild->nalloc);
1085 vpbc = ((interaction_function[ftype].flags & IF_VSITE) &&
1086 vsite->vsite_pbc_loc != NULL);
1089 nral1 = 1 + NRAL(ftype);
1090 ftv = ftype - F_VSITE2;
1091 if ((ild->nr + n)/nral1 > vsite->vsite_pbc_loc_nalloc[ftv])
1093 vsite->vsite_pbc_loc_nalloc[ftv] =
1094 over_alloc_large((ild->nr + n)/nral1);
1095 srenew(vsite->vsite_pbc_loc[ftv],
1096 vsite->vsite_pbc_loc_nalloc[ftv]);
1100 for(s=0; s<nsrc; s++)
1102 ils = &src[s].il[ftype];
1103 for(i=0; i<ils->nr; i++)
1105 ild->iatoms[ild->nr+i] = ils->iatoms[i];
1109 for(i=0; i<ils->nr; i+=nral1)
1111 vsite->vsite_pbc_loc[ftv][(ild->nr+i)/nral1] =
1112 vsite_pbc_t[s][ftv][i/nral1];
1121 /* Position restraints need an additional treatment */
1122 if (dest->il[F_POSRES].nr > 0)
1124 n = dest->il[F_POSRES].nr/2;
1125 if (n > dest->iparams_posres_nalloc)
1127 dest->iparams_posres_nalloc = over_alloc_large(n);
1128 srenew(dest->iparams_posres,dest->iparams_posres_nalloc);
1130 /* Set n to the number of original position restraints in dest */
1131 for(s=0; s<nsrc; s++)
1133 n -= src[s].il[F_POSRES].nr/2;
1135 for(s=0; s<nsrc; s++)
1137 for(i=0; i<src[s].il[F_POSRES].nr/2; i++)
1139 /* Correct the index into iparams_posres */
1140 dest->il[F_POSRES].iatoms[n*2] = n;
1141 /* Copy the position restraint force parameters */
1142 dest->iparams_posres[n] = src[s].iparams_posres[i];
1149 /* This function looks up and assigns bonded interactions for zone iz.
1150 * With thread parallelizing each thread acts on a different atom range:
1151 * at_start to at_end.
1153 static int make_bondeds_zone(gmx_domdec_t *dd,
1154 const gmx_domdec_zones_t *zones,
1155 const gmx_molblock_t *molb,
1156 gmx_bool bRCheckMB,ivec rcheck,gmx_bool bRCheck2B,
1158 int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1159 const t_iparams *ip_in,
1160 t_idef *idef,gmx_vsite_t *vsite,
1162 int *vsite_pbc_nalloc,
1164 int at_start,int at_end)
1166 int i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,d,k;
1168 t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
1169 gmx_bool bBCheck,bUse,bLocal;
1175 const gmx_domdec_ns_ranges_t *izone;
1176 gmx_reverse_top_t *rt;
1179 nizone = zones->nizone;
1180 izone = zones->izone;
1182 rt = dd->reverse_top;
1184 bBCheck = rt->bBCheck;
1190 for(i=at_start; i<at_end; i++)
1192 /* Get the global atom number */
1193 i_gl = dd->gatindex[i];
1194 global_atomnr_to_moltype_ind(rt,i_gl,&mb,&mt,&mol,&i_mol);
1195 /* Check all interactions assigned to this atom */
1196 index = rt->ril_mt[mt].index;
1197 rtil = rt->ril_mt[mt].il;
1199 while (j < index[i_mol+1])
1204 if (ftype == F_SETTLE)
1206 /* Settles are only in the reverse top when they
1207 * operate within a charge group. So we can assign
1208 * them without checks. We do this only for performance
1209 * reasons; it could be handled by the code below.
1213 /* Home zone: add this settle to the local topology */
1214 tiatoms[0] = iatoms[0];
1216 tiatoms[2] = i + iatoms[2] - iatoms[1];
1217 tiatoms[3] = i + iatoms[3] - iatoms[1];
1218 add_ifunc(nral,tiatoms,&idef->il[ftype]);
1223 else if (interaction_function[ftype].flags & IF_VSITE)
1225 /* The vsite construction goes where the vsite itself is */
1228 add_vsite(dd->ga2la,index,rtil,ftype,nral,
1230 iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
1237 tiatoms[0] = iatoms[0];
1241 /* Assign single-body interactions to the home zone */
1246 if (ftype == F_POSRES)
1248 add_posres(mol,i_mol,&molb[mb],tiatoms,ip_in,
1251 else if (ftype == F_FBPOSRES)
1253 add_fbposres(mol,i_mol,&molb[mb],tiatoms,ip_in,
1264 /* This is a two-body interaction, we can assign
1265 * analogous to the non-bonded assignments.
1267 if (!ga2la_get(ga2la,i_gl+iatoms[2]-i_mol,&a_loc,&kz))
1277 /* Check zone interaction assignments */
1278 bUse = ((iz < nizone && iz <= kz &&
1279 izone[iz].j0 <= kz && kz < izone[iz].j1) ||
1280 (kz < nizone && iz > kz &&
1281 izone[kz].j0 <= iz && iz < izone[kz].j1));
1286 /* If necessary check the cgcm distance */
1288 dd_dist2(pbc_null,cg_cm,la2lc,
1289 tiatoms[1],tiatoms[2]) >= rc2)
1298 /* Assign this multi-body bonded interaction to
1299 * the local node if we have all the atoms involved
1300 * (local or communicated) and the minimum zone shift
1301 * in each dimension is zero, for dimensions
1302 * with 2 DD cells an extra check may be necessary.
1307 for(k=1; k<=nral && bUse; k++)
1309 bLocal = ga2la_get(ga2la,i_gl+iatoms[k]-i_mol,
1311 if (!bLocal || kz >= zones->n)
1313 /* We do not have this atom of this interaction
1314 * locally, or it comes from more than one cell
1322 for(d=0; d<DIM; d++)
1324 if (zones->shift[kz][d] == 0)
1336 k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
1339 for(d=0; (d<DIM && bUse); d++)
1341 /* Check if the cg_cm distance falls within
1342 * the cut-off to avoid possible multiple
1343 * assignments of bonded interactions.
1347 dd_dist2(pbc_null,cg_cm,la2lc,
1348 tiatoms[k_zero[d]],tiatoms[k_plus[d]]) >= rc2)
1357 /* Add this interaction to the local topology */
1358 add_ifunc(nral,tiatoms,&idef->il[ftype]);
1359 /* Sum so we can check in global_stat
1360 * if we have everything.
1363 !(interaction_function[ftype].flags & IF_LIMZERO))
1373 return nbonded_local;
1376 static void set_no_exclusions_zone(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1377 int iz,t_blocka *lexcls)
1381 a0 = dd->cgindex[zones->cg_range[iz]];
1382 a1 = dd->cgindex[zones->cg_range[iz+1]];
1384 for(a=a0+1; a<a1+1; a++)
1386 lexcls->index[a] = lexcls->nra;
1390 static int make_exclusions_zone(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1391 const gmx_moltype_t *moltype,
1392 gmx_bool bRCheck,real rc2,
1393 int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1397 int cg_start,int cg_end)
1399 int nizone,n,count,jla0,jla1,jla;
1400 int cg,la0,la1,la,a_gl,mb,mt,mol,a_mol,j,aj_mol;
1401 const t_blocka *excls;
1408 jla0 = dd->cgindex[zones->izone[iz].jcg0];
1409 jla1 = dd->cgindex[zones->izone[iz].jcg1];
1411 /* We set the end index, but note that we might not start at zero here */
1412 lexcls->nr = dd->cgindex[cg_end];
1416 for(cg=cg_start; cg<cg_end; cg++)
1418 /* Here we assume the number of exclusions in one charge group
1419 * is never larger than 1000.
1421 if (n+1000 > lexcls->nalloc_a)
1423 lexcls->nalloc_a = over_alloc_large(n+1000);
1424 srenew(lexcls->a,lexcls->nalloc_a);
1426 la0 = dd->cgindex[cg];
1427 la1 = dd->cgindex[cg+1];
1428 if (GET_CGINFO_EXCL_INTER(cginfo[cg]) ||
1429 !GET_CGINFO_EXCL_INTRA(cginfo[cg]))
1431 /* Copy the exclusions from the global top */
1432 for(la=la0; la<la1; la++) {
1433 lexcls->index[la] = n;
1434 a_gl = dd->gatindex[la];
1435 global_atomnr_to_moltype_ind(dd->reverse_top,a_gl,&mb,&mt,&mol,&a_mol);
1436 excls = &moltype[mt].excls;
1437 for(j=excls->index[a_mol]; j<excls->index[a_mol+1]; j++)
1439 aj_mol = excls->a[j];
1440 /* This computation of jla is only correct intra-cg */
1441 jla = la + aj_mol - a_mol;
1442 if (jla >= la0 && jla < la1)
1444 /* This is an intra-cg exclusion. We can skip
1445 * the global indexing and distance checking.
1447 /* Intra-cg exclusions are only required
1448 * for the home zone.
1452 lexcls->a[n++] = jla;
1453 /* Check to avoid double counts */
1462 /* This is a inter-cg exclusion */
1463 /* Since exclusions are pair interactions,
1464 * just like non-bonded interactions,
1465 * they can be assigned properly up
1466 * to the DD cutoff (not cutoff_min as
1467 * for the other bonded interactions).
1469 if (ga2la_get(ga2la,a_gl+aj_mol-a_mol,&jla,&cell))
1471 if (iz == 0 && cell == 0)
1473 lexcls->a[n++] = jla;
1474 /* Check to avoid double counts */
1480 else if (jla >= jla0 && jla < jla1 &&
1482 dd_dist2(pbc_null,cg_cm,la2lc,la,jla) < rc2))
1484 /* jla > la, since jla0 > la */
1485 lexcls->a[n++] = jla;
1495 /* There are no inter-cg excls and this cg is self-excluded.
1496 * These exclusions are only required for zone 0,
1497 * since other zones do not see themselves.
1501 for(la=la0; la<la1; la++)
1503 lexcls->index[la] = n;
1504 for(j=la0; j<la1; j++)
1509 count += ((la1 - la0)*(la1 - la0 - 1))/2;
1513 /* We don't need exclusions for this cg */
1514 for(la=la0; la<la1; la++)
1516 lexcls->index[la] = n;
1522 lexcls->index[lexcls->nr] = n;
1528 static void check_alloc_index(t_blocka *ba,int nindex_max)
1530 if (nindex_max+1 > ba->nalloc_index)
1532 ba->nalloc_index = over_alloc_dd(nindex_max+1);
1533 srenew(ba->index,ba->nalloc_index);
1537 static void check_exclusions_alloc(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1543 nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1545 check_alloc_index(lexcls,nr);
1547 for(thread=1; thread<dd->reverse_top->nthread; thread++)
1549 check_alloc_index(&dd->reverse_top->excl_thread[thread],nr);
1553 static void finish_local_exclusions(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1558 lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1560 if (dd->n_intercg_excl == 0)
1562 /* There are no exclusions involving non-home charge groups,
1563 * but we need to set the indices for neighborsearching.
1565 la0 = dd->cgindex[zones->izone[0].cg1];
1566 for(la=la0; la<lexcls->nr; la++)
1568 lexcls->index[la] = lexcls->nra;
1571 /* nr is only used to loop over the exclusions for Ewald and RF,
1572 * so we can set it to the number of home atoms for efficiency.
1574 lexcls->nr = dd->cgindex[zones->izone[0].cg1];
1578 static void clear_idef(t_idef *idef)
1582 /* Clear the counts */
1583 for(ftype=0; ftype<F_NRE; ftype++)
1585 idef->il[ftype].nr = 0;
1589 static int make_local_bondeds_excls(gmx_domdec_t *dd,
1590 gmx_domdec_zones_t *zones,
1591 const gmx_mtop_t *mtop,
1593 gmx_bool bRCheckMB,ivec rcheck,gmx_bool bRCheck2B,
1595 int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1596 t_idef *idef,gmx_vsite_t *vsite,
1597 t_blocka *lexcls,int *excl_count)
1599 int nzone_bondeds,nzone_excl;
1604 gmx_reverse_top_t *rt;
1606 if (dd->reverse_top->bMultiCGmols)
1608 nzone_bondeds = zones->n;
1612 /* Only single charge group molecules, so interactions don't
1613 * cross zone boundaries and we only need to assign in the home zone.
1618 if (dd->n_intercg_excl > 0)
1620 /* We only use exclusions from i-zones to i- and j-zones */
1621 nzone_excl = zones->nizone;
1625 /* There are no inter-cg exclusions and only zone 0 sees itself */
1629 check_exclusions_alloc(dd,zones,lexcls);
1631 rt = dd->reverse_top;
1635 /* Clear the counts */
1643 for(iz=0; iz<nzone_bondeds; iz++)
1645 cg0 = zones->cg_range[iz];
1646 cg1 = zones->cg_range[iz+1];
1648 #pragma omp parallel for num_threads(rt->nthread) schedule(static)
1649 for(thread=0; thread<rt->nthread; thread++)
1655 int *vsite_pbc_nalloc;
1658 cg0t = cg0 + ((cg1 - cg0)* thread )/rt->nthread;
1659 cg1t = cg0 + ((cg1 - cg0)*(thread+1))/rt->nthread;
1667 idef_t = &rt->idef_thread[thread];
1671 if (vsite && vsite->n_intercg_vsite > 0)
1675 vsite_pbc = vsite->vsite_pbc_loc;
1676 vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
1680 vsite_pbc = rt->vsite_pbc[thread];
1681 vsite_pbc_nalloc = rt->vsite_pbc_nalloc[thread];
1687 vsite_pbc_nalloc = NULL;
1690 rt->nbonded_thread[thread] =
1691 make_bondeds_zone(dd,zones,
1693 bRCheckMB,rcheck,bRCheck2B,rc2,
1694 la2lc,pbc_null,cg_cm,idef->iparams,
1696 vsite,vsite_pbc,vsite_pbc_nalloc,
1698 dd->cgindex[cg0t],dd->cgindex[cg1t]);
1700 if (iz < nzone_excl)
1708 excl_t = &rt->excl_thread[thread];
1713 rt->excl_count_thread[thread] =
1714 make_exclusions_zone(dd,zones,
1715 mtop->moltype,bRCheck2B,rc2,
1716 la2lc,pbc_null,cg_cm,cginfo,
1723 if (rt->nthread > 1)
1725 combine_idef(idef,rt->idef_thread+1,rt->nthread-1,
1726 vsite,rt->vsite_pbc+1);
1729 for(thread=0; thread<rt->nthread; thread++)
1731 nbonded_local += rt->nbonded_thread[thread];
1734 if (iz < nzone_excl)
1736 if (rt->nthread > 1)
1738 combine_blocka(lexcls,rt->excl_thread+1,rt->nthread-1);
1741 for(thread=0; thread<rt->nthread; thread++)
1743 *excl_count += rt->excl_count_thread[thread];
1748 /* Some zones might not have exclusions, but some code still needs to
1749 * loop over the index, so we set the indices here.
1751 for(iz=nzone_excl; iz<zones->nizone; iz++)
1753 set_no_exclusions_zone(dd,zones,iz,lexcls);
1756 finish_local_exclusions(dd,zones,lexcls);
1759 fprintf(debug,"We have %d exclusions, check count %d\n",
1760 lexcls->nra,*excl_count);
1763 return nbonded_local;
1766 void dd_make_local_cgs(gmx_domdec_t *dd,t_block *lcgs)
1768 lcgs->nr = dd->ncg_tot;
1769 lcgs->index = dd->cgindex;
1772 void dd_make_local_top(FILE *fplog,
1773 gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1774 int npbcdim,matrix box,
1775 rvec cellsize_min,ivec npulse,
1779 gmx_mtop_t *mtop,gmx_localtop_t *ltop)
1781 gmx_bool bUniqueExcl,bRCheckMB,bRCheck2B,bRCheckExcl;
1785 t_pbc pbc,*pbc_null=NULL;
1789 fprintf(debug,"Making local topology\n");
1792 dd_make_local_cgs(dd,<op->cgs);
1796 bRCheckExcl = FALSE;
1798 if (dd->reverse_top->bMultiCGmols)
1800 /* We need to check to which cell bondeds should be assigned */
1801 rc = dd_cutoff_twobody(dd);
1804 fprintf(debug,"Two-body bonded cut-off distance is %g\n",rc);
1807 /* Should we check cg_cm distances when assigning bonded interactions? */
1808 for(d=0; d<DIM; d++)
1811 /* Only need to check for dimensions where the part of the box
1812 * that is not communicated is smaller than the cut-off.
1814 if (d < npbcdim && dd->nc[d] > 1 &&
1815 (dd->nc[d] - npulse[d])*cellsize_min[d] < 2*rc)
1822 /* Check for interactions between two atoms,
1823 * where we can allow interactions up to the cut-off,
1824 * instead of up to the smallest cell dimension.
1831 "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
1832 d,cellsize_min[d],d,rcheck[d],bRCheck2B);
1835 if (dd->reverse_top->bExclRequired)
1837 bRCheckExcl = bRCheck2B;
1841 /* If we don't have forces on exclusions,
1842 * we don't care about exclusions being assigned mulitple times.
1844 bRCheckExcl = FALSE;
1846 if (bRCheckMB || bRCheck2B)
1851 set_pbc_dd(&pbc,fr->ePBC,dd,TRUE,box);
1862 make_local_bondeds_excls(dd,zones,mtop,fr->cginfo,
1863 bRCheckMB,rcheck,bRCheck2B,rc,
1867 <op->excls,&nexcl);
1869 /* The ilist is not sorted yet,
1870 * we can only do this when we have the charge arrays.
1872 ltop->idef.ilsort = ilsortUNKNOWN;
1874 if (dd->reverse_top->bExclRequired)
1876 dd->nbonded_local += nexcl;
1878 forcerec_set_excl_load(fr,ltop,NULL);
1881 ltop->atomtypes = mtop->atomtypes;
1883 /* For an error message only */
1884 dd->reverse_top->err_top_global = mtop;
1885 dd->reverse_top->err_top_local = ltop;
1888 void dd_sort_local_top(gmx_domdec_t *dd,t_mdatoms *mdatoms,
1889 gmx_localtop_t *ltop)
1891 if (dd->reverse_top->ilsort == ilsortNO_FE)
1893 ltop->idef.ilsort = ilsortNO_FE;
1897 gmx_sort_ilist_fe(<op->idef,mdatoms->chargeA,mdatoms->chargeB);
1901 gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global)
1903 gmx_localtop_t *top;
1908 top->idef.ntypes = top_global->ffparams.ntypes;
1909 top->idef.atnr = top_global->ffparams.atnr;
1910 top->idef.functype = top_global->ffparams.functype;
1911 top->idef.iparams = top_global->ffparams.iparams;
1912 top->idef.fudgeQQ = top_global->ffparams.fudgeQQ;
1913 top->idef.cmap_grid= top_global->ffparams.cmap_grid;
1915 for(i=0; i<F_NRE; i++)
1917 top->idef.il[i].iatoms = NULL;
1918 top->idef.il[i].nalloc = 0;
1920 top->idef.ilsort = ilsortUNKNOWN;
1925 void dd_init_local_state(gmx_domdec_t *dd,
1926 t_state *state_global,t_state *state_local)
1928 int buf[NITEM_DD_INIT_LOCAL_STATE];
1932 buf[0] = state_global->flags;
1933 buf[1] = state_global->ngtc;
1934 buf[2] = state_global->nnhpres;
1935 buf[3] = state_global->nhchainlength;
1936 buf[4] = state_global->dfhist.nlambda;
1938 dd_bcast(dd,NITEM_DD_INIT_LOCAL_STATE*sizeof(int),buf);
1940 init_state(state_local,0,buf[1],buf[2],buf[3],buf[4]);
1941 state_local->flags = buf[0];
1943 /* With Langevin Dynamics we need to make proper storage space
1944 * in the global and local state for the random numbers.
1946 if (state_local->flags & (1<<estLD_RNG))
1948 if (DDMASTER(dd) && state_global->nrngi > 1)
1950 state_global->nrng = dd->nnodes*gmx_rng_n();
1951 srenew(state_global->ld_rng,state_global->nrng);
1953 state_local->nrng = gmx_rng_n();
1954 snew(state_local->ld_rng,state_local->nrng);
1956 if (state_local->flags & (1<<estLD_RNGI))
1958 if (DDMASTER(dd) && state_global->nrngi > 1)
1960 state_global->nrngi = dd->nnodes;
1961 srenew(state_global->ld_rngi,state_global->nrngi);
1963 snew(state_local->ld_rngi,1);
1967 static void check_link(t_blocka *link,int cg_gl,int cg_gl_j)
1973 for(k=link->index[cg_gl]; k<link->index[cg_gl+1]; k++)
1975 if (link->a[k] == cg_gl_j)
1982 /* Add this charge group link */
1983 if (link->index[cg_gl+1]+1 > link->nalloc_a)
1985 link->nalloc_a = over_alloc_large(link->index[cg_gl+1]+1);
1986 srenew(link->a,link->nalloc_a);
1988 link->a[link->index[cg_gl+1]] = cg_gl_j;
1989 link->index[cg_gl+1]++;
1993 static int *make_at2cg(t_block *cgs)
1997 snew(at2cg,cgs->index[cgs->nr]);
1998 for(cg=0; cg<cgs->nr; cg++)
2000 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
2009 t_blocka *make_charge_group_links(gmx_mtop_t *mtop,gmx_domdec_t *dd,
2010 cginfo_mb_t *cginfo_mb)
2012 gmx_reverse_top_t *rt;
2013 int mb,cg_offset,cg,cg_gl,a,aj,i,j,ftype,nral,nlink_mol,mol,ncgi;
2014 gmx_molblock_t *molb;
2015 gmx_moltype_t *molt;
2019 gmx_reverse_ilist_t ril;
2021 cginfo_mb_t *cgi_mb;
2023 /* For each charge group make a list of other charge groups
2024 * in the system that a linked to it via bonded interactions
2025 * which are also stored in reverse_top.
2028 rt = dd->reverse_top;
2031 snew(link->index,ncg_mtop(mtop)+1);
2038 for(mb=0; mb<mtop->nmolblock; mb++)
2040 molb = &mtop->molblock[mb];
2041 if (molb->nmol == 0)
2045 molt = &mtop->moltype[molb->type];
2047 excls = &molt->excls;
2048 a2c = make_at2cg(cgs);
2049 /* Make a reverse ilist in which the interactions are linked
2050 * to all atoms, not only the first atom as in gmx_reverse_top.
2051 * The constraints are discarded here.
2053 make_reverse_ilist(molt,NULL,FALSE,FALSE,FALSE,TRUE,&ril);
2055 cgi_mb = &cginfo_mb[mb];
2057 for(cg=0; cg<cgs->nr; cg++)
2059 cg_gl = cg_offset + cg;
2060 link->index[cg_gl+1] = link->index[cg_gl];
2061 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
2064 while (i < ril.index[a+1])
2066 ftype = ril.il[i++];
2068 /* Skip the ifunc index */
2070 for(j=0; j<nral; j++)
2075 check_link(link,cg_gl,cg_offset+a2c[aj]);
2078 i += nral_rt(ftype);
2080 if (rt->bExclRequired)
2082 /* Exclusions always go both ways */
2083 for(j=excls->index[a]; j<excls->index[a+1]; j++)
2088 check_link(link,cg_gl,cg_offset+a2c[aj]);
2093 if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
2095 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
2099 nlink_mol = link->index[cg_offset+cgs->nr] - link->index[cg_offset];
2101 cg_offset += cgs->nr;
2103 destroy_reverse_ilist(&ril);
2108 fprintf(debug,"molecule type '%s' %d cgs has %d cg links through bonded interac.\n",*molt->name,cgs->nr,nlink_mol);
2113 /* Copy the data for the rest of the molecules in this block */
2114 link->nalloc_a += (molb->nmol - 1)*nlink_mol;
2115 srenew(link->a,link->nalloc_a);
2116 for(mol=1; mol<molb->nmol; mol++)
2118 for(cg=0; cg<cgs->nr; cg++)
2120 cg_gl = cg_offset + cg;
2121 link->index[cg_gl+1] =
2122 link->index[cg_gl+1-cgs->nr] + nlink_mol;
2123 for(j=link->index[cg_gl]; j<link->index[cg_gl+1]; j++)
2125 link->a[j] = link->a[j-nlink_mol] + cgs->nr;
2127 if (link->index[cg_gl+1] - link->index[cg_gl] > 0 &&
2128 cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
2130 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
2134 cg_offset += cgs->nr;
2141 fprintf(debug,"Of the %d charge groups %d are linked via bonded interactions\n",ncg_mtop(mtop),ncgi);
2147 static void bonded_cg_distance_mol(gmx_moltype_t *molt,int *at2cg,
2148 gmx_bool bBCheck,gmx_bool bExcl,rvec *cg_cm,
2149 real *r_2b,int *ft2b,int *a2_1,int *a2_2,
2150 real *r_mb,int *ftmb,int *am_1,int *am_2)
2152 int ftype,nral,i,j,ai,aj,cgi,cgj;
2155 real r2_2b,r2_mb,rij2;
2159 for(ftype=0; ftype<F_NRE; ftype++)
2161 if (dd_check_ftype(ftype,bBCheck,FALSE,FALSE))
2163 il = &molt->ilist[ftype];
2167 for(i=0; i<il->nr; i+=1+nral)
2169 for(ai=0; ai<nral; ai++)
2171 cgi = at2cg[il->iatoms[i+1+ai]];
2172 for(aj=0; aj<nral; aj++) {
2173 cgj = at2cg[il->iatoms[i+1+aj]];
2176 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
2177 if (nral == 2 && rij2 > r2_2b)
2181 *a2_1 = il->iatoms[i+1+ai];
2182 *a2_2 = il->iatoms[i+1+aj];
2184 if (nral > 2 && rij2 > r2_mb)
2188 *am_1 = il->iatoms[i+1+ai];
2189 *am_2 = il->iatoms[i+1+aj];
2200 excls = &molt->excls;
2201 for(ai=0; ai<excls->nr; ai++)
2204 for(j=excls->index[ai]; j<excls->index[ai+1]; j++) {
2205 cgj = at2cg[excls->a[j]];
2208 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
2218 *r_2b = sqrt(r2_2b);
2219 *r_mb = sqrt(r2_mb);
2222 static void get_cgcm_mol(gmx_moltype_t *molt,gmx_ffparams_t *ffparams,
2223 int ePBC,t_graph *graph,matrix box,
2225 rvec *x,rvec *xs,rvec *cg_cm)
2229 if (ePBC != epbcNONE)
2231 mk_mshift(NULL,graph,ePBC,box,x);
2233 shift_x(graph,box,x,xs);
2234 /* By doing an extra mk_mshift the molecules that are broken
2235 * because they were e.g. imported from another software
2236 * will be made whole again. Such are the healing powers
2239 mk_mshift(NULL,graph,ePBC,box,xs);
2243 /* We copy the coordinates so the original coordinates remain
2244 * unchanged, just to be 100% sure that we do not affect
2245 * binary reproducibility of simulations.
2247 n = molt->cgs.index[molt->cgs.nr];
2250 copy_rvec(x[i],xs[i]);
2256 construct_vsites(NULL,vsite,xs,NULL,0.0,NULL,
2257 ffparams->iparams,molt->ilist,
2258 epbcNONE,TRUE,NULL,NULL,NULL);
2261 calc_cgcm(NULL,0,molt->cgs.nr,&molt->cgs,xs,cg_cm);
2264 static int have_vsite_molt(gmx_moltype_t *molt)
2270 for(i=0; i<F_NRE; i++)
2272 if ((interaction_function[i].flags & IF_VSITE) &&
2273 molt->ilist[i].nr > 0) {
2281 void dd_bonded_cg_distance(FILE *fplog,
2282 gmx_domdec_t *dd,gmx_mtop_t *mtop,
2283 t_inputrec *ir,rvec *x,matrix box,
2285 real *r_2b,real *r_mb)
2287 gmx_bool bExclRequired;
2288 int mb,cg_offset,at_offset,*at2cg,mol;
2291 gmx_molblock_t *molb;
2292 gmx_moltype_t *molt;
2294 real rmol_2b,rmol_mb;
2295 int ft2b=-1,a_2b_1=-1,a_2b_2=-1,ftmb=-1,a_mb_1=-1,a_mb_2=-1;
2296 int ftm2b=-1,amol_2b_1=-1,amol_2b_2=-1,ftmmb=-1,amol_mb_1=-1,amol_mb_2=-1;
2298 bExclRequired = IR_EXCL_FORCES(*ir);
2300 /* For gmx_vsite_t everything 0 should work (without pbc) */
2307 for(mb=0; mb<mtop->nmolblock; mb++)
2309 molb = &mtop->molblock[mb];
2310 molt = &mtop->moltype[molb->type];
2311 if (molt->cgs.nr == 1 || molb->nmol == 0)
2313 cg_offset += molb->nmol*molt->cgs.nr;
2314 at_offset += molb->nmol*molt->atoms.nr;
2318 if (ir->ePBC != epbcNONE)
2320 mk_graph_ilist(NULL,molt->ilist,0,molt->atoms.nr,FALSE,FALSE,
2324 at2cg = make_at2cg(&molt->cgs);
2325 snew(xs,molt->atoms.nr);
2326 snew(cg_cm,molt->cgs.nr);
2327 for(mol=0; mol<molb->nmol; mol++)
2329 get_cgcm_mol(molt,&mtop->ffparams,ir->ePBC,&graph,box,
2330 have_vsite_molt(molt) ? vsite : NULL,
2331 x+at_offset,xs,cg_cm);
2333 bonded_cg_distance_mol(molt,at2cg,bBCheck,bExclRequired,cg_cm,
2334 &rmol_2b,&ftm2b,&amol_2b_1,&amol_2b_2,
2335 &rmol_mb,&ftmmb,&amol_mb_1,&amol_mb_2);
2336 if (rmol_2b > *r_2b)
2340 a_2b_1 = at_offset + amol_2b_1;
2341 a_2b_2 = at_offset + amol_2b_2;
2343 if (rmol_mb > *r_mb)
2347 a_mb_1 = at_offset + amol_mb_1;
2348 a_mb_2 = at_offset + amol_mb_2;
2351 cg_offset += molt->cgs.nr;
2352 at_offset += molt->atoms.nr;
2357 if (ir->ePBC != epbcNONE)
2366 if (fplog && (ft2b >= 0 || ftmb >= 0))
2369 "Initial maximum inter charge-group distances:\n");
2373 " two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2374 *r_2b,interaction_function[ft2b].longname,
2380 " multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2381 *r_mb,interaction_function[ftmb].longname,