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_vsite(gmx_ga2la_t ga2la
,int *index
,int *rtil
,
832 gmx_bool bHomeA
,int a
,int a_gl
,int a_mol
,
834 t_idef
*idef
,int **vsite_pbc
,int *vsite_pbc_nalloc
)
836 int k
,ak_gl
,vsi
,pbc_a_mol
;
837 t_iatom tiatoms
[1+MAXATOMLIST
],*iatoms_r
;
838 int j
,ftype_r
,nral_r
;
841 tiatoms
[0] = iatoms
[0];
845 /* We know the local index of the first atom */
850 /* Convert later in make_local_vsites */
851 tiatoms
[1] = -a_gl
- 1;
854 for(k
=2; k
<1+nral
; k
++)
856 ak_gl
= a_gl
+ iatoms
[k
] - a_mol
;
857 if (!ga2la_get_home(ga2la
,ak_gl
,&tiatoms
[k
]))
859 /* Copy the global index, convert later in make_local_vsites */
860 tiatoms
[k
] = -(ak_gl
+ 1);
864 /* Add this interaction to the local topology */
865 add_ifunc(nral
,tiatoms
,&idef
->il
[ftype
]);
868 vsi
= idef
->il
[ftype
].nr
/(1+nral
) - 1;
869 if (vsi
>= vsite_pbc_nalloc
[ftype
-F_VSITE2
])
871 vsite_pbc_nalloc
[ftype
-F_VSITE2
] = over_alloc_large(vsi
+1);
872 srenew(vsite_pbc
[ftype
-F_VSITE2
],vsite_pbc_nalloc
[ftype
-F_VSITE2
]);
876 pbc_a_mol
= iatoms
[1+nral
+1];
879 /* The pbc flag is one of the following two options:
880 * -2: vsite and all constructing atoms are within the same cg, no pbc
881 * -1: vsite and its first constructing atom are in the same cg, do pbc
883 vsite_pbc
[ftype
-F_VSITE2
][vsi
] = pbc_a_mol
;
887 /* Set the pbc atom for this vsite so we can make its pbc
888 * identical to the rest of the atoms in its charge group.
889 * Since the order of the atoms does not change within a charge
890 * group, we do not need the global to local atom index.
892 vsite_pbc
[ftype
-F_VSITE2
][vsi
] = a
+ pbc_a_mol
- iatoms
[1];
897 /* This vsite is non-home (required for recursion),
898 * and therefore there is no charge group to match pbc with.
899 * But we always turn on full_pbc to assure that higher order
900 * recursion works correctly.
902 vsite_pbc
[ftype
-F_VSITE2
][vsi
] = -1;
908 /* Check for recursion */
909 for(k
=2; k
<1+nral
; k
++)
911 if ((iatoms
[1+nral
] & (2<<k
)) && (tiatoms
[k
] < 0))
913 /* This construction atoms is a vsite and not a home atom */
916 fprintf(debug
,"Constructing atom %d of vsite atom %d is a vsite and non-home\n",iatoms
[k
]+1,a_mol
+1);
918 /* Find the vsite construction */
920 /* Check all interactions assigned to this atom */
921 j
= index
[iatoms
[k
]];
922 while (j
< index
[iatoms
[k
]+1])
925 nral_r
= NRAL(ftype_r
);
926 if (interaction_function
[ftype_r
].flags
& IF_VSITE
)
928 /* Add this vsite (recursion) */
929 add_vsite(ga2la
,index
,rtil
,ftype_r
,nral_r
,
930 FALSE
,-1,a_gl
+iatoms
[k
]-iatoms
[1],iatoms
[k
],
931 rtil
+j
,idef
,vsite_pbc
,vsite_pbc_nalloc
);
944 static void make_la2lc(gmx_domdec_t
*dd
)
946 int *cgindex
,*la2lc
,cg
,a
;
948 cgindex
= dd
->cgindex
;
950 if (dd
->nat_tot
> dd
->la2lc_nalloc
)
952 dd
->la2lc_nalloc
= over_alloc_dd(dd
->nat_tot
);
953 snew(dd
->la2lc
,dd
->la2lc_nalloc
);
957 /* Make the local atom to local cg index */
958 for(cg
=0; cg
<dd
->ncg_tot
; cg
++)
960 for(a
=cgindex
[cg
]; a
<cgindex
[cg
+1]; a
++)
967 static real
dd_dist2(t_pbc
*pbc_null
,rvec
*cg_cm
,const int *la2lc
,int i
,int j
)
973 pbc_dx_aiuc(pbc_null
,cg_cm
[la2lc
[i
]],cg_cm
[la2lc
[j
]],dx
);
977 rvec_sub(cg_cm
[la2lc
[i
]],cg_cm
[la2lc
[j
]],dx
);
983 /* Append the nsrc t_blocka block structures in src to *dest */
984 static void combine_blocka(t_blocka
*dest
,const t_blocka
*src
,int nsrc
)
990 for(s
=0; s
<nsrc
; s
++)
994 if (ni
+ 1 > dest
->nalloc_index
)
996 dest
->nalloc_index
= over_alloc_large(ni
+1);
997 srenew(dest
->index
,dest
->nalloc_index
);
999 if (dest
->nra
+ na
> dest
->nalloc_a
)
1001 dest
->nalloc_a
= over_alloc_large(dest
->nra
+na
);
1002 srenew(dest
->a
,dest
->nalloc_a
);
1004 for(s
=0; s
<nsrc
; s
++)
1006 for(i
=dest
->nr
+1; i
<src
[s
].nr
+1; i
++)
1008 dest
->index
[i
] = dest
->nra
+ src
[s
].index
[i
];
1010 for(i
=0; i
<src
[s
].nra
; i
++)
1012 dest
->a
[dest
->nra
+i
] = src
[s
].a
[i
];
1014 dest
->nr
= src
[s
].nr
;
1015 dest
->nra
+= src
[s
].nra
;
1019 /* Append the nsrc t_idef structures in src to *dest,
1020 * virtual sites need special attention, as pbc info differs per vsite.
1022 static void combine_idef(t_idef
*dest
,const t_idef
*src
,int nsrc
,
1023 gmx_vsite_t
*vsite
,int ***vsite_pbc_t
)
1031 for(ftype
=0; ftype
<F_NRE
; ftype
++)
1034 for(s
=0; s
<nsrc
; s
++)
1036 n
+= src
[s
].il
[ftype
].nr
;
1040 ild
= &dest
->il
[ftype
];
1042 if (ild
->nr
+ n
> ild
->nalloc
)
1044 ild
->nalloc
= over_alloc_large(ild
->nr
+n
);
1045 srenew(ild
->iatoms
,ild
->nalloc
);
1048 vpbc
= ((interaction_function
[ftype
].flags
& IF_VSITE
) &&
1049 vsite
->vsite_pbc_loc
!= NULL
);
1052 nral1
= 1 + NRAL(ftype
);
1053 ftv
= ftype
- F_VSITE2
;
1054 if ((ild
->nr
+ n
)/nral1
> vsite
->vsite_pbc_loc_nalloc
[ftv
])
1056 vsite
->vsite_pbc_loc_nalloc
[ftv
] =
1057 over_alloc_large((ild
->nr
+ n
)/nral1
);
1058 srenew(vsite
->vsite_pbc_loc
[ftv
],
1059 vsite
->vsite_pbc_loc_nalloc
[ftv
]);
1063 for(s
=0; s
<nsrc
; s
++)
1065 ils
= &src
[s
].il
[ftype
];
1066 for(i
=0; i
<ils
->nr
; i
++)
1068 ild
->iatoms
[ild
->nr
+i
] = ils
->iatoms
[i
];
1072 for(i
=0; i
<ils
->nr
; i
+=nral1
)
1074 vsite
->vsite_pbc_loc
[ftv
][(ild
->nr
+i
)/nral1
] =
1075 vsite_pbc_t
[s
][ftv
][i
/nral1
];
1084 /* Position restraints need an additional treatment */
1085 if (dest
->il
[F_POSRES
].nr
> 0)
1087 n
= dest
->il
[F_POSRES
].nr
/2;
1088 if (n
> dest
->iparams_posres_nalloc
)
1090 dest
->iparams_posres_nalloc
= over_alloc_large(n
);
1091 srenew(dest
->iparams_posres
,dest
->iparams_posres_nalloc
);
1093 /* Set n to the number of original position restraints in dest */
1094 for(s
=0; s
<nsrc
; s
++)
1096 n
-= src
[s
].il
[F_POSRES
].nr
/2;
1098 for(s
=0; s
<nsrc
; s
++)
1100 for(i
=0; i
<src
[s
].il
[F_POSRES
].nr
/2; i
++)
1102 /* Correct the index into iparams_posres */
1103 dest
->il
[F_POSRES
].iatoms
[n
*2] = n
;
1104 /* Copy the position restraint force parameters */
1105 dest
->iparams_posres
[n
] = src
[s
].iparams_posres
[i
];
1112 /* This function looks up and assigns bonded interactions for zone iz.
1113 * With thread parallelizing each thread acts on a different atom range:
1114 * at_start to at_end.
1116 static int make_bondeds_zone(gmx_domdec_t
*dd
,
1117 const gmx_domdec_zones_t
*zones
,
1118 const gmx_molblock_t
*molb
,
1119 gmx_bool bRCheckMB
,ivec rcheck
,gmx_bool bRCheck2B
,
1121 int *la2lc
,t_pbc
*pbc_null
,rvec
*cg_cm
,
1122 const t_iparams
*ip_in
,
1123 t_idef
*idef
,gmx_vsite_t
*vsite
,
1125 int *vsite_pbc_nalloc
,
1127 int at_start
,int at_end
)
1129 int i
,i_gl
,mb
,mt
,mol
,i_mol
,j
,ftype
,nral
,d
,k
;
1131 t_iatom
*iatoms
,tiatoms
[1+MAXATOMLIST
];
1132 gmx_bool bBCheck
,bUse
,bLocal
;
1138 const gmx_domdec_ns_ranges_t
*izone
;
1139 gmx_reverse_top_t
*rt
;
1142 nizone
= zones
->nizone
;
1143 izone
= zones
->izone
;
1145 rt
= dd
->reverse_top
;
1147 bBCheck
= rt
->bBCheck
;
1153 for(i
=at_start
; i
<at_end
; i
++)
1155 /* Get the global atom number */
1156 i_gl
= dd
->gatindex
[i
];
1157 global_atomnr_to_moltype_ind(rt
,i_gl
,&mb
,&mt
,&mol
,&i_mol
);
1158 /* Check all interactions assigned to this atom */
1159 index
= rt
->ril_mt
[mt
].index
;
1160 rtil
= rt
->ril_mt
[mt
].il
;
1162 while (j
< index
[i_mol
+1])
1167 if (ftype
== F_SETTLE
)
1169 /* Settles are only in the reverse top when they
1170 * operate within a charge group. So we can assign
1171 * them without checks. We do this only for performance
1172 * reasons; it could be handled by the code below.
1176 /* Home zone: add this settle to the local topology */
1177 tiatoms
[0] = iatoms
[0];
1179 tiatoms
[2] = i
+ iatoms
[2] - iatoms
[1];
1180 tiatoms
[3] = i
+ iatoms
[3] - iatoms
[1];
1181 add_ifunc(nral
,tiatoms
,&idef
->il
[ftype
]);
1186 else if (interaction_function
[ftype
].flags
& IF_VSITE
)
1188 /* The vsite construction goes where the vsite itself is */
1191 add_vsite(dd
->ga2la
,index
,rtil
,ftype
,nral
,
1193 iatoms
,idef
,vsite_pbc
,vsite_pbc_nalloc
);
1200 tiatoms
[0] = iatoms
[0];
1204 /* Assign single-body interactions to the home zone */
1209 if (ftype
== F_POSRES
)
1211 add_posres(mol
,i_mol
,&molb
[mb
],tiatoms
,ip_in
,
1222 /* This is a two-body interaction, we can assign
1223 * analogous to the non-bonded assignments.
1225 if (!ga2la_get(ga2la
,i_gl
+iatoms
[2]-i_mol
,&a_loc
,&kz
))
1235 /* Check zone interaction assignments */
1236 bUse
= ((iz
< nizone
&& iz
<= kz
&&
1237 izone
[iz
].j0
<= kz
&& kz
< izone
[iz
].j1
) ||
1238 (kz
< nizone
&& iz
> kz
&&
1239 izone
[kz
].j0
<= iz
&& iz
< izone
[kz
].j1
));
1244 /* If necessary check the cgcm distance */
1246 dd_dist2(pbc_null
,cg_cm
,la2lc
,
1247 tiatoms
[1],tiatoms
[2]) >= rc2
)
1256 /* Assign this multi-body bonded interaction to
1257 * the local node if we have all the atoms involved
1258 * (local or communicated) and the minimum zone shift
1259 * in each dimension is zero, for dimensions
1260 * with 2 DD cells an extra check may be necessary.
1265 for(k
=1; k
<=nral
&& bUse
; k
++)
1267 bLocal
= ga2la_get(ga2la
,i_gl
+iatoms
[k
]-i_mol
,
1269 if (!bLocal
|| kz
>= zones
->n
)
1271 /* We do not have this atom of this interaction
1272 * locally, or it comes from more than one cell
1280 for(d
=0; d
<DIM
; d
++)
1282 if (zones
->shift
[kz
][d
] == 0)
1294 k_zero
[XX
] && k_zero
[YY
] && k_zero
[ZZ
]);
1297 for(d
=0; (d
<DIM
&& bUse
); d
++)
1299 /* Check if the cg_cm distance falls within
1300 * the cut-off to avoid possible multiple
1301 * assignments of bonded interactions.
1305 dd_dist2(pbc_null
,cg_cm
,la2lc
,
1306 tiatoms
[k_zero
[d
]],tiatoms
[k_plus
[d
]]) >= rc2
)
1315 /* Add this interaction to the local topology */
1316 add_ifunc(nral
,tiatoms
,&idef
->il
[ftype
]);
1317 /* Sum so we can check in global_stat
1318 * if we have everything.
1321 !(interaction_function
[ftype
].flags
& IF_LIMZERO
))
1331 return nbonded_local
;
1334 static void set_no_exclusions_zone(gmx_domdec_t
*dd
,gmx_domdec_zones_t
*zones
,
1335 int iz
,t_blocka
*lexcls
)
1339 a0
= dd
->cgindex
[zones
->cg_range
[iz
]];
1340 a1
= dd
->cgindex
[zones
->cg_range
[iz
+1]];
1342 for(a
=a0
+1; a
<a1
+1; a
++)
1344 lexcls
->index
[a
] = lexcls
->nra
;
1348 static int make_exclusions_zone(gmx_domdec_t
*dd
,gmx_domdec_zones_t
*zones
,
1349 const gmx_moltype_t
*moltype
,
1350 gmx_bool bRCheck
,real rc2
,
1351 int *la2lc
,t_pbc
*pbc_null
,rvec
*cg_cm
,
1355 int cg_start
,int cg_end
)
1357 int nizone
,n
,count
,jla0
,jla1
,jla
;
1358 int cg
,la0
,la1
,la
,a_gl
,mb
,mt
,mol
,a_mol
,j
,aj_mol
;
1359 const t_blocka
*excls
;
1366 jla0
= dd
->cgindex
[zones
->izone
[iz
].jcg0
];
1367 jla1
= dd
->cgindex
[zones
->izone
[iz
].jcg1
];
1369 /* We set the end index, but note that we might not start at zero here */
1370 lexcls
->nr
= dd
->cgindex
[cg_end
];
1374 for(cg
=cg_start
; cg
<cg_end
; cg
++)
1376 /* Here we assume the number of exclusions in one charge group
1377 * is never larger than 1000.
1379 if (n
+1000 > lexcls
->nalloc_a
)
1381 lexcls
->nalloc_a
= over_alloc_large(n
+1000);
1382 srenew(lexcls
->a
,lexcls
->nalloc_a
);
1384 la0
= dd
->cgindex
[cg
];
1385 la1
= dd
->cgindex
[cg
+1];
1386 if (GET_CGINFO_EXCL_INTER(cginfo
[cg
]) ||
1387 !GET_CGINFO_EXCL_INTRA(cginfo
[cg
]))
1389 /* Copy the exclusions from the global top */
1390 for(la
=la0
; la
<la1
; la
++) {
1391 lexcls
->index
[la
] = n
;
1392 a_gl
= dd
->gatindex
[la
];
1393 global_atomnr_to_moltype_ind(dd
->reverse_top
,a_gl
,&mb
,&mt
,&mol
,&a_mol
);
1394 excls
= &moltype
[mt
].excls
;
1395 for(j
=excls
->index
[a_mol
]; j
<excls
->index
[a_mol
+1]; j
++)
1397 aj_mol
= excls
->a
[j
];
1398 /* This computation of jla is only correct intra-cg */
1399 jla
= la
+ aj_mol
- a_mol
;
1400 if (jla
>= la0
&& jla
< la1
)
1402 /* This is an intra-cg exclusion. We can skip
1403 * the global indexing and distance checking.
1405 /* Intra-cg exclusions are only required
1406 * for the home zone.
1410 lexcls
->a
[n
++] = jla
;
1411 /* Check to avoid double counts */
1420 /* This is a inter-cg exclusion */
1421 /* Since exclusions are pair interactions,
1422 * just like non-bonded interactions,
1423 * they can be assigned properly up
1424 * to the DD cutoff (not cutoff_min as
1425 * for the other bonded interactions).
1427 if (ga2la_get(ga2la
,a_gl
+aj_mol
-a_mol
,&jla
,&cell
))
1429 if (iz
== 0 && cell
== 0)
1431 lexcls
->a
[n
++] = jla
;
1432 /* Check to avoid double counts */
1438 else if (jla
>= jla0
&& jla
< jla1
&&
1440 dd_dist2(pbc_null
,cg_cm
,la2lc
,la
,jla
) < rc2
))
1442 /* jla > la, since jla0 > la */
1443 lexcls
->a
[n
++] = jla
;
1453 /* There are no inter-cg excls and this cg is self-excluded.
1454 * These exclusions are only required for zone 0,
1455 * since other zones do not see themselves.
1459 for(la
=la0
; la
<la1
; la
++)
1461 lexcls
->index
[la
] = n
;
1462 for(j
=la0
; j
<la1
; j
++)
1467 count
+= ((la1
- la0
)*(la1
- la0
- 1))/2;
1471 /* We don't need exclusions for this cg */
1472 for(la
=la0
; la
<la1
; la
++)
1474 lexcls
->index
[la
] = n
;
1480 lexcls
->index
[lexcls
->nr
] = n
;
1486 static void check_alloc_index(t_blocka
*ba
,int nindex_max
)
1488 if (nindex_max
+1 > ba
->nalloc_index
)
1490 ba
->nalloc_index
= over_alloc_dd(nindex_max
+1);
1491 srenew(ba
->index
,ba
->nalloc_index
);
1495 static void check_exclusions_alloc(gmx_domdec_t
*dd
,gmx_domdec_zones_t
*zones
,
1501 nr
= dd
->cgindex
[zones
->izone
[zones
->nizone
-1].cg1
];
1503 check_alloc_index(lexcls
,nr
);
1505 for(thread
=1; thread
<dd
->reverse_top
->nthread
; thread
++)
1507 check_alloc_index(&dd
->reverse_top
->excl_thread
[thread
],nr
);
1511 static void finish_local_exclusions(gmx_domdec_t
*dd
,gmx_domdec_zones_t
*zones
,
1516 lexcls
->nr
= dd
->cgindex
[zones
->izone
[zones
->nizone
-1].cg1
];
1518 if (dd
->n_intercg_excl
== 0)
1520 /* There are no exclusions involving non-home charge groups,
1521 * but we need to set the indices for neighborsearching.
1523 la0
= dd
->cgindex
[zones
->izone
[0].cg1
];
1524 for(la
=la0
; la
<lexcls
->nr
; la
++)
1526 lexcls
->index
[la
] = lexcls
->nra
;
1529 /* nr is only used to loop over the exclusions for Ewald and RF,
1530 * so we can set it to the number of home atoms for efficiency.
1532 lexcls
->nr
= dd
->cgindex
[zones
->izone
[0].cg1
];
1536 static void clear_idef(t_idef
*idef
)
1540 /* Clear the counts */
1541 for(ftype
=0; ftype
<F_NRE
; ftype
++)
1543 idef
->il
[ftype
].nr
= 0;
1547 static int make_local_bondeds_excls(gmx_domdec_t
*dd
,
1548 gmx_domdec_zones_t
*zones
,
1549 const gmx_mtop_t
*mtop
,
1551 gmx_bool bRCheckMB
,ivec rcheck
,gmx_bool bRCheck2B
,
1553 int *la2lc
,t_pbc
*pbc_null
,rvec
*cg_cm
,
1554 t_idef
*idef
,gmx_vsite_t
*vsite
,
1555 t_blocka
*lexcls
,int *excl_count
)
1557 int nzone_bondeds
,nzone_excl
;
1562 gmx_reverse_top_t
*rt
;
1564 if (dd
->reverse_top
->bMultiCGmols
)
1566 nzone_bondeds
= zones
->n
;
1570 /* Only single charge group molecules, so interactions don't
1571 * cross zone boundaries and we only need to assign in the home zone.
1576 if (dd
->n_intercg_excl
> 0)
1578 /* We only use exclusions from i-zones to i- and j-zones */
1579 nzone_excl
= zones
->nizone
;
1583 /* There are no inter-cg exclusions and only zone 0 sees itself */
1587 check_exclusions_alloc(dd
,zones
,lexcls
);
1589 rt
= dd
->reverse_top
;
1593 /* Clear the counts */
1601 for(iz
=0; iz
<nzone_bondeds
; iz
++)
1603 cg0
= zones
->cg_range
[iz
];
1604 cg1
= zones
->cg_range
[iz
+1];
1606 #pragma omp parallel for num_threads(rt->nthread) schedule(static)
1607 for(thread
=0; thread
<rt
->nthread
; thread
++)
1613 int *vsite_pbc_nalloc
;
1616 cg0t
= cg0
+ ((cg1
- cg0
)* thread
)/rt
->nthread
;
1617 cg1t
= cg0
+ ((cg1
- cg0
)*(thread
+1))/rt
->nthread
;
1625 idef_t
= &rt
->idef_thread
[thread
];
1629 if (vsite
&& vsite
->n_intercg_vsite
> 0)
1633 vsite_pbc
= vsite
->vsite_pbc_loc
;
1634 vsite_pbc_nalloc
= vsite
->vsite_pbc_loc_nalloc
;
1638 vsite_pbc
= rt
->vsite_pbc
[thread
];
1639 vsite_pbc_nalloc
= rt
->vsite_pbc_nalloc
[thread
];
1645 vsite_pbc_nalloc
= NULL
;
1648 rt
->nbonded_thread
[thread
] =
1649 make_bondeds_zone(dd
,zones
,
1651 bRCheckMB
,rcheck
,bRCheck2B
,rc2
,
1652 la2lc
,pbc_null
,cg_cm
,idef
->iparams
,
1654 vsite
,vsite_pbc
,vsite_pbc_nalloc
,
1656 dd
->cgindex
[cg0t
],dd
->cgindex
[cg1t
]);
1658 if (iz
< nzone_excl
)
1666 excl_t
= &rt
->excl_thread
[thread
];
1671 rt
->excl_count_thread
[thread
] =
1672 make_exclusions_zone(dd
,zones
,
1673 mtop
->moltype
,bRCheck2B
,rc2
,
1674 la2lc
,pbc_null
,cg_cm
,cginfo
,
1681 if (rt
->nthread
> 1)
1683 combine_idef(idef
,rt
->idef_thread
+1,rt
->nthread
-1,
1684 vsite
,rt
->vsite_pbc
+1);
1687 for(thread
=0; thread
<rt
->nthread
; thread
++)
1689 nbonded_local
+= rt
->nbonded_thread
[thread
];
1692 if (iz
< nzone_excl
)
1694 if (rt
->nthread
> 1)
1696 combine_blocka(lexcls
,rt
->excl_thread
+1,rt
->nthread
-1);
1699 for(thread
=0; thread
<rt
->nthread
; thread
++)
1701 *excl_count
+= rt
->excl_count_thread
[thread
];
1706 /* Some zones might not have exclusions, but some code still needs to
1707 * loop over the index, so we set the indices here.
1709 for(iz
=nzone_excl
; iz
<zones
->nizone
; iz
++)
1711 set_no_exclusions_zone(dd
,zones
,iz
,lexcls
);
1714 finish_local_exclusions(dd
,zones
,lexcls
);
1717 fprintf(debug
,"We have %d exclusions, check count %d\n",
1718 lexcls
->nra
,*excl_count
);
1721 return nbonded_local
;
1724 void dd_make_local_cgs(gmx_domdec_t
*dd
,t_block
*lcgs
)
1726 lcgs
->nr
= dd
->ncg_tot
;
1727 lcgs
->index
= dd
->cgindex
;
1730 void dd_make_local_top(FILE *fplog
,
1731 gmx_domdec_t
*dd
,gmx_domdec_zones_t
*zones
,
1732 int npbcdim
,matrix box
,
1733 rvec cellsize_min
,ivec npulse
,
1737 gmx_mtop_t
*mtop
,gmx_localtop_t
*ltop
)
1739 gmx_bool bUniqueExcl
,bRCheckMB
,bRCheck2B
,bRCheckExcl
;
1743 t_pbc pbc
,*pbc_null
=NULL
;
1747 fprintf(debug
,"Making local topology\n");
1750 dd_make_local_cgs(dd
,<op
->cgs
);
1754 bRCheckExcl
= FALSE
;
1756 if (dd
->reverse_top
->bMultiCGmols
)
1758 /* We need to check to which cell bondeds should be assigned */
1759 rc
= dd_cutoff_twobody(dd
);
1762 fprintf(debug
,"Two-body bonded cut-off distance is %g\n",rc
);
1765 /* Should we check cg_cm distances when assigning bonded interactions? */
1766 for(d
=0; d
<DIM
; d
++)
1769 /* Only need to check for dimensions where the part of the box
1770 * that is not communicated is smaller than the cut-off.
1772 if (d
< npbcdim
&& dd
->nc
[d
] > 1 &&
1773 (dd
->nc
[d
] - npulse
[d
])*cellsize_min
[d
] < 2*rc
)
1780 /* Check for interactions between two atoms,
1781 * where we can allow interactions up to the cut-off,
1782 * instead of up to the smallest cell dimension.
1789 "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
1790 d
,cellsize_min
[d
],d
,rcheck
[d
],bRCheck2B
);
1793 if (dd
->reverse_top
->bExclRequired
)
1795 bRCheckExcl
= bRCheck2B
;
1799 /* If we don't have forces on exclusions,
1800 * we don't care about exclusions being assigned mulitple times.
1802 bRCheckExcl
= FALSE
;
1804 if (bRCheckMB
|| bRCheck2B
)
1809 set_pbc_dd(&pbc
,fr
->ePBC
,dd
,TRUE
,box
);
1820 make_local_bondeds_excls(dd
,zones
,mtop
,fr
->cginfo
,
1821 bRCheckMB
,rcheck
,bRCheck2B
,rc
,
1825 <op
->excls
,&nexcl
);
1827 /* The ilist is not sorted yet,
1828 * we can only do this when we have the charge arrays.
1830 ltop
->idef
.ilsort
= ilsortUNKNOWN
;
1832 if (dd
->reverse_top
->bExclRequired
)
1834 dd
->nbonded_local
+= nexcl
;
1836 forcerec_set_excl_load(fr
,ltop
,NULL
);
1839 ltop
->atomtypes
= mtop
->atomtypes
;
1841 /* For an error message only */
1842 dd
->reverse_top
->err_top_global
= mtop
;
1843 dd
->reverse_top
->err_top_local
= ltop
;
1846 void dd_sort_local_top(gmx_domdec_t
*dd
,t_mdatoms
*mdatoms
,
1847 gmx_localtop_t
*ltop
)
1849 if (dd
->reverse_top
->ilsort
== ilsortNO_FE
)
1851 ltop
->idef
.ilsort
= ilsortNO_FE
;
1855 gmx_sort_ilist_fe(<op
->idef
,mdatoms
->chargeA
,mdatoms
->chargeB
);
1859 gmx_localtop_t
*dd_init_local_top(gmx_mtop_t
*top_global
)
1861 gmx_localtop_t
*top
;
1866 top
->idef
.ntypes
= top_global
->ffparams
.ntypes
;
1867 top
->idef
.atnr
= top_global
->ffparams
.atnr
;
1868 top
->idef
.functype
= top_global
->ffparams
.functype
;
1869 top
->idef
.iparams
= top_global
->ffparams
.iparams
;
1870 top
->idef
.fudgeQQ
= top_global
->ffparams
.fudgeQQ
;
1871 top
->idef
.cmap_grid
= top_global
->ffparams
.cmap_grid
;
1873 for(i
=0; i
<F_NRE
; i
++)
1875 top
->idef
.il
[i
].iatoms
= NULL
;
1876 top
->idef
.il
[i
].nalloc
= 0;
1878 top
->idef
.ilsort
= ilsortUNKNOWN
;
1883 void dd_init_local_state(gmx_domdec_t
*dd
,
1884 t_state
*state_global
,t_state
*state_local
)
1886 int buf
[NITEM_DD_INIT_LOCAL_STATE
];
1890 buf
[0] = state_global
->flags
;
1891 buf
[1] = state_global
->ngtc
;
1892 buf
[2] = state_global
->nnhpres
;
1893 buf
[3] = state_global
->nhchainlength
;
1894 buf
[4] = state_global
->dfhist
.nlambda
;
1896 dd_bcast(dd
,NITEM_DD_INIT_LOCAL_STATE
*sizeof(int),buf
);
1898 init_state(state_local
,0,buf
[1],buf
[2],buf
[3],buf
[4]);
1899 state_local
->flags
= buf
[0];
1901 /* With Langevin Dynamics we need to make proper storage space
1902 * in the global and local state for the random numbers.
1904 if (state_local
->flags
& (1<<estLD_RNG
))
1906 if (DDMASTER(dd
) && state_global
->nrngi
> 1)
1908 state_global
->nrng
= dd
->nnodes
*gmx_rng_n();
1909 srenew(state_global
->ld_rng
,state_global
->nrng
);
1911 state_local
->nrng
= gmx_rng_n();
1912 snew(state_local
->ld_rng
,state_local
->nrng
);
1914 if (state_local
->flags
& (1<<estLD_RNGI
))
1916 if (DDMASTER(dd
) && state_global
->nrngi
> 1)
1918 state_global
->nrngi
= dd
->nnodes
;
1919 srenew(state_global
->ld_rngi
,state_global
->nrngi
);
1921 snew(state_local
->ld_rngi
,1);
1925 static void check_link(t_blocka
*link
,int cg_gl
,int cg_gl_j
)
1931 for(k
=link
->index
[cg_gl
]; k
<link
->index
[cg_gl
+1]; k
++)
1933 if (link
->a
[k
] == cg_gl_j
)
1940 /* Add this charge group link */
1941 if (link
->index
[cg_gl
+1]+1 > link
->nalloc_a
)
1943 link
->nalloc_a
= over_alloc_large(link
->index
[cg_gl
+1]+1);
1944 srenew(link
->a
,link
->nalloc_a
);
1946 link
->a
[link
->index
[cg_gl
+1]] = cg_gl_j
;
1947 link
->index
[cg_gl
+1]++;
1951 static int *make_at2cg(t_block
*cgs
)
1955 snew(at2cg
,cgs
->index
[cgs
->nr
]);
1956 for(cg
=0; cg
<cgs
->nr
; cg
++)
1958 for(a
=cgs
->index
[cg
]; a
<cgs
->index
[cg
+1]; a
++)
1967 t_blocka
*make_charge_group_links(gmx_mtop_t
*mtop
,gmx_domdec_t
*dd
,
1968 cginfo_mb_t
*cginfo_mb
)
1970 gmx_reverse_top_t
*rt
;
1971 int mb
,cg_offset
,cg
,cg_gl
,a
,aj
,i
,j
,ftype
,nral
,nlink_mol
,mol
,ncgi
;
1972 gmx_molblock_t
*molb
;
1973 gmx_moltype_t
*molt
;
1977 gmx_reverse_ilist_t ril
;
1979 cginfo_mb_t
*cgi_mb
;
1981 /* For each charge group make a list of other charge groups
1982 * in the system that a linked to it via bonded interactions
1983 * which are also stored in reverse_top.
1986 rt
= dd
->reverse_top
;
1989 snew(link
->index
,ncg_mtop(mtop
)+1);
1996 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
1998 molb
= &mtop
->molblock
[mb
];
1999 if (molb
->nmol
== 0)
2003 molt
= &mtop
->moltype
[molb
->type
];
2005 excls
= &molt
->excls
;
2006 a2c
= make_at2cg(cgs
);
2007 /* Make a reverse ilist in which the interactions are linked
2008 * to all atoms, not only the first atom as in gmx_reverse_top.
2009 * The constraints are discarded here.
2011 make_reverse_ilist(molt
,NULL
,FALSE
,FALSE
,FALSE
,TRUE
,&ril
);
2013 cgi_mb
= &cginfo_mb
[mb
];
2015 for(cg
=0; cg
<cgs
->nr
; cg
++)
2017 cg_gl
= cg_offset
+ cg
;
2018 link
->index
[cg_gl
+1] = link
->index
[cg_gl
];
2019 for(a
=cgs
->index
[cg
]; a
<cgs
->index
[cg
+1]; a
++)
2022 while (i
< ril
.index
[a
+1])
2024 ftype
= ril
.il
[i
++];
2026 /* Skip the ifunc index */
2028 for(j
=0; j
<nral
; j
++)
2033 check_link(link
,cg_gl
,cg_offset
+a2c
[aj
]);
2036 i
+= nral_rt(ftype
);
2038 if (rt
->bExclRequired
)
2040 /* Exclusions always go both ways */
2041 for(j
=excls
->index
[a
]; j
<excls
->index
[a
+1]; j
++)
2046 check_link(link
,cg_gl
,cg_offset
+a2c
[aj
]);
2051 if (link
->index
[cg_gl
+1] - link
->index
[cg_gl
] > 0)
2053 SET_CGINFO_BOND_INTER(cgi_mb
->cginfo
[cg
]);
2057 nlink_mol
= link
->index
[cg_offset
+cgs
->nr
] - link
->index
[cg_offset
];
2059 cg_offset
+= cgs
->nr
;
2061 destroy_reverse_ilist(&ril
);
2066 fprintf(debug
,"molecule type '%s' %d cgs has %d cg links through bonded interac.\n",*molt
->name
,cgs
->nr
,nlink_mol
);
2071 /* Copy the data for the rest of the molecules in this block */
2072 link
->nalloc_a
+= (molb
->nmol
- 1)*nlink_mol
;
2073 srenew(link
->a
,link
->nalloc_a
);
2074 for(mol
=1; mol
<molb
->nmol
; mol
++)
2076 for(cg
=0; cg
<cgs
->nr
; cg
++)
2078 cg_gl
= cg_offset
+ cg
;
2079 link
->index
[cg_gl
+1] =
2080 link
->index
[cg_gl
+1-cgs
->nr
] + nlink_mol
;
2081 for(j
=link
->index
[cg_gl
]; j
<link
->index
[cg_gl
+1]; j
++)
2083 link
->a
[j
] = link
->a
[j
-nlink_mol
] + cgs
->nr
;
2085 if (link
->index
[cg_gl
+1] - link
->index
[cg_gl
] > 0 &&
2086 cg_gl
- cgi_mb
->cg_start
< cgi_mb
->cg_mod
)
2088 SET_CGINFO_BOND_INTER(cgi_mb
->cginfo
[cg_gl
- cgi_mb
->cg_start
]);
2092 cg_offset
+= cgs
->nr
;
2099 fprintf(debug
,"Of the %d charge groups %d are linked via bonded interactions\n",ncg_mtop(mtop
),ncgi
);
2105 static void bonded_cg_distance_mol(gmx_moltype_t
*molt
,int *at2cg
,
2106 gmx_bool bBCheck
,gmx_bool bExcl
,rvec
*cg_cm
,
2107 real
*r_2b
,int *ft2b
,int *a2_1
,int *a2_2
,
2108 real
*r_mb
,int *ftmb
,int *am_1
,int *am_2
)
2110 int ftype
,nral
,i
,j
,ai
,aj
,cgi
,cgj
;
2113 real r2_2b
,r2_mb
,rij2
;
2117 for(ftype
=0; ftype
<F_NRE
; ftype
++)
2119 if (dd_check_ftype(ftype
,bBCheck
,FALSE
,FALSE
))
2121 il
= &molt
->ilist
[ftype
];
2125 for(i
=0; i
<il
->nr
; i
+=1+nral
)
2127 for(ai
=0; ai
<nral
; ai
++)
2129 cgi
= at2cg
[il
->iatoms
[i
+1+ai
]];
2130 for(aj
=0; aj
<nral
; aj
++) {
2131 cgj
= at2cg
[il
->iatoms
[i
+1+aj
]];
2134 rij2
= distance2(cg_cm
[cgi
],cg_cm
[cgj
]);
2135 if (nral
== 2 && rij2
> r2_2b
)
2139 *a2_1
= il
->iatoms
[i
+1+ai
];
2140 *a2_2
= il
->iatoms
[i
+1+aj
];
2142 if (nral
> 2 && rij2
> r2_mb
)
2146 *am_1
= il
->iatoms
[i
+1+ai
];
2147 *am_2
= il
->iatoms
[i
+1+aj
];
2158 excls
= &molt
->excls
;
2159 for(ai
=0; ai
<excls
->nr
; ai
++)
2162 for(j
=excls
->index
[ai
]; j
<excls
->index
[ai
+1]; j
++) {
2163 cgj
= at2cg
[excls
->a
[j
]];
2166 rij2
= distance2(cg_cm
[cgi
],cg_cm
[cgj
]);
2176 *r_2b
= sqrt(r2_2b
);
2177 *r_mb
= sqrt(r2_mb
);
2180 static void get_cgcm_mol(gmx_moltype_t
*molt
,gmx_ffparams_t
*ffparams
,
2181 int ePBC
,t_graph
*graph
,matrix box
,
2183 rvec
*x
,rvec
*xs
,rvec
*cg_cm
)
2187 if (ePBC
!= epbcNONE
)
2189 mk_mshift(NULL
,graph
,ePBC
,box
,x
);
2191 shift_x(graph
,box
,x
,xs
);
2192 /* By doing an extra mk_mshift the molecules that are broken
2193 * because they were e.g. imported from another software
2194 * will be made whole again. Such are the healing powers
2197 mk_mshift(NULL
,graph
,ePBC
,box
,xs
);
2201 /* We copy the coordinates so the original coordinates remain
2202 * unchanged, just to be 100% sure that we do not affect
2203 * binary reproducibility of simulations.
2205 n
= molt
->cgs
.index
[molt
->cgs
.nr
];
2208 copy_rvec(x
[i
],xs
[i
]);
2214 construct_vsites(NULL
,vsite
,xs
,NULL
,0.0,NULL
,
2215 ffparams
->iparams
,molt
->ilist
,
2216 epbcNONE
,TRUE
,NULL
,NULL
,NULL
);
2219 calc_cgcm(NULL
,0,molt
->cgs
.nr
,&molt
->cgs
,xs
,cg_cm
);
2222 static int have_vsite_molt(gmx_moltype_t
*molt
)
2228 for(i
=0; i
<F_NRE
; i
++)
2230 if ((interaction_function
[i
].flags
& IF_VSITE
) &&
2231 molt
->ilist
[i
].nr
> 0) {
2239 void dd_bonded_cg_distance(FILE *fplog
,
2240 gmx_domdec_t
*dd
,gmx_mtop_t
*mtop
,
2241 t_inputrec
*ir
,rvec
*x
,matrix box
,
2243 real
*r_2b
,real
*r_mb
)
2245 gmx_bool bExclRequired
;
2246 int mb
,cg_offset
,at_offset
,*at2cg
,mol
;
2249 gmx_molblock_t
*molb
;
2250 gmx_moltype_t
*molt
;
2252 real rmol_2b
,rmol_mb
;
2253 int ft2b
=-1,a_2b_1
=-1,a_2b_2
=-1,ftmb
=-1,a_mb_1
=-1,a_mb_2
=-1;
2254 int ftm2b
=-1,amol_2b_1
=-1,amol_2b_2
=-1,ftmmb
=-1,amol_mb_1
=-1,amol_mb_2
=-1;
2256 bExclRequired
= IR_EXCL_FORCES(*ir
);
2258 /* For gmx_vsite_t everything 0 should work (without pbc) */
2265 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
2267 molb
= &mtop
->molblock
[mb
];
2268 molt
= &mtop
->moltype
[molb
->type
];
2269 if (molt
->cgs
.nr
== 1 || molb
->nmol
== 0)
2271 cg_offset
+= molb
->nmol
*molt
->cgs
.nr
;
2272 at_offset
+= molb
->nmol
*molt
->atoms
.nr
;
2276 if (ir
->ePBC
!= epbcNONE
)
2278 mk_graph_ilist(NULL
,molt
->ilist
,0,molt
->atoms
.nr
,FALSE
,FALSE
,
2282 at2cg
= make_at2cg(&molt
->cgs
);
2283 snew(xs
,molt
->atoms
.nr
);
2284 snew(cg_cm
,molt
->cgs
.nr
);
2285 for(mol
=0; mol
<molb
->nmol
; mol
++)
2287 get_cgcm_mol(molt
,&mtop
->ffparams
,ir
->ePBC
,&graph
,box
,
2288 have_vsite_molt(molt
) ? vsite
: NULL
,
2289 x
+at_offset
,xs
,cg_cm
);
2291 bonded_cg_distance_mol(molt
,at2cg
,bBCheck
,bExclRequired
,cg_cm
,
2292 &rmol_2b
,&ftm2b
,&amol_2b_1
,&amol_2b_2
,
2293 &rmol_mb
,&ftmmb
,&amol_mb_1
,&amol_mb_2
);
2294 if (rmol_2b
> *r_2b
)
2298 a_2b_1
= at_offset
+ amol_2b_1
;
2299 a_2b_2
= at_offset
+ amol_2b_2
;
2301 if (rmol_mb
> *r_mb
)
2305 a_mb_1
= at_offset
+ amol_mb_1
;
2306 a_mb_2
= at_offset
+ amol_mb_2
;
2309 cg_offset
+= molt
->cgs
.nr
;
2310 at_offset
+= molt
->atoms
.nr
;
2315 if (ir
->ePBC
!= epbcNONE
)
2324 if (fplog
&& (ft2b
>= 0 || ftmb
>= 0))
2327 "Initial maximum inter charge-group distances:\n");
2331 " two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2332 *r_2b
,interaction_function
[ft2b
].longname
,
2338 " multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2339 *r_mb
,interaction_function
[ftmb
].longname
,