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
28 #include "domdec_network.h"
29 #include "mtop_util.h"
30 #include "gmx_ga2la.h"
39 typedef struct gmx_domdec_specat_comm
{
40 /* The atom indices we need from the surrounding cells */
44 /* The number of indices to receive during the setup */
46 /* The atoms to send */
47 gmx_specatsend_t spas
[DIM
][2];
57 /* The range in the local buffer(s) for received atoms */
60 } gmx_domdec_specat_comm_t
;
62 typedef struct gmx_domdec_constraints
{
65 /* The fully local and connected constraints */
67 /* The global constraint number, only required for clearing gc_req */
71 /* Boolean that tells if a global constraint index has been requested */
73 /* Global to local communicated constraint atom only index */
75 } gmx_domdec_constraints_t
;
78 static void dd_move_f_specat(gmx_domdec_t
*dd
,gmx_domdec_specat_comm_t
*spac
,
81 gmx_specatsend_t
*spas
;
83 int n
,n0
,n1
,d
,dim
,dir
,i
;
89 for(d
=dd
->ndim
-1; d
>=0; d
--)
94 /* Pulse the grid forward and backward */
100 /* Send and receive the coordinates */
101 dd_sendrecv2_rvec(dd
,d
,
102 f
+n
+n1
,n0
,vbuf
,spas
[0].nsend
,
103 f
+n
,n1
,vbuf
+spas
[0].nsend
,spas
[1].nsend
);
104 for(dir
=0; dir
<2; dir
++)
106 bPBC
= ((dir
== 0 && dd
->ci
[dim
] == 0) ||
107 (dir
== 1 && dd
->ci
[dim
] == dd
->nc
[dim
]-1));
108 bScrew
= (bPBC
&& dd
->bScrewPBC
&& dim
== XX
);
110 spas
= &spac
->spas
[d
][dir
];
111 /* Sum the buffer into the required forces */
112 if (!bPBC
|| (!bScrew
&& fshift
== NULL
))
114 for(i
=0; i
<spas
->nsend
; i
++)
116 rvec_inc(f
[spas
->a
[i
]],*vbuf
);
123 vis
[dim
] = (dir
==0 ? 1 : -1);
127 /* Sum and add to shift forces */
128 for(i
=0; i
<spas
->nsend
; i
++)
130 rvec_inc(f
[spas
->a
[i
]],*vbuf
);
131 rvec_inc(fshift
[is
],*vbuf
);
137 /* Rotate the forces */
138 for(i
=0; i
<spas
->nsend
; i
++)
140 f
[spas
->a
[i
]][XX
] += (*vbuf
)[XX
];
141 f
[spas
->a
[i
]][YY
] -= (*vbuf
)[YY
];
142 f
[spas
->a
[i
]][ZZ
] -= (*vbuf
)[ZZ
];
145 rvec_inc(fshift
[is
],*vbuf
);
155 /* Two cells, so we only need to communicate one way */
156 spas
= &spac
->spas
[d
][0];
158 /* Send and receive the coordinates */
159 dd_sendrecv_rvec(dd
,d
,dddirForward
,
160 f
+n
,spas
->nrecv
,spac
->vbuf
,spas
->nsend
);
161 /* Sum the buffer into the required forces */
162 if (dd
->bScrewPBC
&& dim
== XX
&&
164 dd
->ci
[dim
] == dd
->nc
[dim
]-1))
166 for(i
=0; i
<spas
->nsend
; i
++)
168 /* Rotate the force */
169 f
[spas
->a
[i
]][XX
] += spac
->vbuf
[i
][XX
];
170 f
[spas
->a
[i
]][YY
] -= spac
->vbuf
[i
][YY
];
171 f
[spas
->a
[i
]][ZZ
] -= spac
->vbuf
[i
][ZZ
];
176 for(i
=0; i
<spas
->nsend
; i
++)
178 rvec_inc(f
[spas
->a
[i
]],spac
->vbuf
[i
]);
185 void dd_move_f_vsites(gmx_domdec_t
*dd
,rvec
*f
,rvec
*fshift
)
189 dd_move_f_specat(dd
,dd
->vsite_comm
,f
,fshift
);
193 void dd_clear_f_vsites(gmx_domdec_t
*dd
,rvec
*f
)
199 for(i
=dd
->vsite_comm
->at_start
; i
<dd
->vsite_comm
->at_end
; i
++)
206 static void dd_move_x_specat(gmx_domdec_t
*dd
,gmx_domdec_specat_comm_t
*spac
,
207 matrix box
,rvec
*x0
,rvec
*x1
)
209 gmx_specatsend_t
*spas
;
211 int nvec
,v
,n
,nn
,ns0
,ns1
,nr0
,nr1
,nr
,d
,dim
,dir
,i
;
212 gmx_bool bPBC
,bScrew
=FALSE
;
222 for(d
=0; d
<dd
->ndim
; d
++)
227 /* Pulse the grid forward and backward */
229 for(dir
=0; dir
<2; dir
++)
231 if (dir
== 0 && dd
->ci
[dim
] == 0)
234 bScrew
= (dd
->bScrewPBC
&& dim
== XX
);
235 copy_rvec(box
[dim
],shift
);
237 else if (dir
== 1 && dd
->ci
[dim
] == dd
->nc
[dim
]-1)
240 bScrew
= (dd
->bScrewPBC
&& dim
== XX
);
243 shift
[i
] = -box
[dim
][i
];
251 spas
= &spac
->spas
[d
][dir
];
252 for(v
=0; v
<nvec
; v
++)
254 x
= (v
== 0 ? x0
: x1
);
255 /* Copy the required coordinates to the send buffer */
259 for(i
=0; i
<spas
->nsend
; i
++)
261 copy_rvec(x
[spas
->a
[i
]],*vbuf
);
267 /* Shift coordinates */
268 for(i
=0; i
<spas
->nsend
; i
++)
270 rvec_add(x
[spas
->a
[i
]],shift
,*vbuf
);
276 /* Shift and rotate coordinates */
277 for(i
=0; i
<spas
->nsend
; i
++)
279 (*vbuf
)[XX
] = x
[spas
->a
[i
]][XX
] + shift
[XX
];
280 (*vbuf
)[YY
] = box
[YY
][YY
] - x
[spas
->a
[i
]][YY
] + shift
[YY
];
281 (*vbuf
)[ZZ
] = box
[ZZ
][ZZ
] - x
[spas
->a
[i
]][ZZ
] + shift
[ZZ
];
287 /* Send and receive the coordinates */
288 spas
= spac
->spas
[d
];
295 dd_sendrecv2_rvec(dd
,d
,
296 spac
->vbuf
+ns0
,ns1
,x0
+n
,nr1
,
297 spac
->vbuf
,ns0
,x0
+n
+nr1
,nr0
);
301 /* Communicate both vectors in one buffer */
303 dd_sendrecv2_rvec(dd
,d
,
304 spac
->vbuf
+2*ns0
,2*ns1
,rbuf
,2*nr1
,
305 spac
->vbuf
,2*ns0
,rbuf
+2*nr1
,2*nr0
);
306 /* Split the buffer into the two vectors */
308 for(dir
=1; dir
>=0; dir
--)
310 nr
= spas
[dir
].nrecv
;
313 x
= (v
== 0 ? x0
: x1
);
316 copy_rvec(*rbuf
,x
[nn
+i
]);
327 spas
= &spac
->spas
[d
][0];
328 /* Copy the required coordinates to the send buffer */
330 for(v
=0; v
<nvec
; v
++)
332 x
= (v
== 0 ? x0
: x1
);
333 if (dd
->bScrewPBC
&& dim
== XX
&&
334 (dd
->ci
[XX
] == 0 || dd
->ci
[XX
] == dd
->nc
[XX
]-1))
336 /* Here we only perform the rotation, the rest of the pbc
337 * is handled in the constraint or viste routines.
339 for(i
=0; i
<spas
->nsend
; i
++)
341 (*vbuf
)[XX
] = x
[spas
->a
[i
]][XX
];
342 (*vbuf
)[YY
] = box
[YY
][YY
] - x
[spas
->a
[i
]][YY
];
343 (*vbuf
)[ZZ
] = box
[ZZ
][ZZ
] - x
[spas
->a
[i
]][ZZ
];
349 for(i
=0; i
<spas
->nsend
; i
++)
351 copy_rvec(x
[spas
->a
[i
]],*vbuf
);
356 /* Send and receive the coordinates */
359 dd_sendrecv_rvec(dd
,d
,dddirBackward
,
360 spac
->vbuf
,spas
->nsend
,x0
+n
,spas
->nrecv
);
364 /* Communicate both vectors in one buffer */
366 dd_sendrecv_rvec(dd
,d
,dddirBackward
,
367 spac
->vbuf
,2*spas
->nsend
,rbuf
,2*spas
->nrecv
);
368 /* Split the buffer into the two vectors */
372 x
= (v
== 0 ? x0
: x1
);
375 copy_rvec(*rbuf
,x
[n
+i
]);
385 void dd_move_x_constraints(gmx_domdec_t
*dd
,matrix box
,rvec
*x0
,rvec
*x1
)
387 if (dd
->constraint_comm
)
389 dd_move_x_specat(dd
,dd
->constraint_comm
,box
,x0
,x1
);
393 void dd_move_x_vsites(gmx_domdec_t
*dd
,matrix box
,rvec
*x
)
397 dd_move_x_specat(dd
,dd
->vsite_comm
,box
,x
,NULL
);
401 int *dd_constraints_nlocalatoms(gmx_domdec_t
*dd
)
405 return dd
->constraints
->con_nlocat
;
413 void dd_clear_local_constraint_indices(gmx_domdec_t
*dd
)
415 gmx_domdec_constraints_t
*dc
;
418 dc
= dd
->constraints
;
420 for(i
=0; i
<dc
->ncon
; i
++)
422 dc
->gc_req
[dc
->con_gl
[i
]] = 0;
425 if (dd
->constraint_comm
)
427 for(i
=dd
->constraint_comm
->at_start
; i
<dd
->constraint_comm
->at_end
; i
++)
429 dc
->ga2la
[dd
->gatindex
[i
]] = -1;
434 void dd_clear_local_vsite_indices(gmx_domdec_t
*dd
)
440 for(i
=dd
->vsite_comm
->at_start
; i
<dd
->vsite_comm
->at_end
; i
++)
442 dd
->ga2la_vsite
[dd
->gatindex
[i
]] = -1;
447 static int setup_specat_communication(gmx_domdec_t
*dd
,
448 gmx_domdec_specat_comm_t
*spac
,
452 const char *specat_type
,
455 int nsend
[2],nlast
,nsend_zero
[2]={0,0},*nsend_ptr
;
456 int d
,dim
,ndir
,dir
,nr
,ns
,i
,nrecv_local
,n0
,start
,ireq
,ind
,buf
[2];
457 int nat_tot_specat
,nat_tot_prev
,nalloc_old
;
458 gmx_bool bPBC
,bFirst
;
459 gmx_specatsend_t
*spas
;
463 fprintf(debug
,"Begin setup_specat_communication for %s\n",specat_type
);
466 /* nsend[0]: the number of atoms requested by this node only,
467 * we communicate this for more efficients checks
468 * nsend[1]: the total number of requested atoms
470 nsend
[0] = spac
->nind_req
;
473 for(d
=dd
->ndim
-1; d
>=0; d
--)
475 /* Pulse the grid forward and backward */
477 bPBC
= (dim
< dd
->npbcdim
);
478 if (dd
->nc
[dim
] == 2)
480 /* Only 2 cells, so we only need to communicate once */
487 for(dir
=0; dir
<ndir
; dir
++)
491 ((dir
== 0 && dd
->ci
[dim
] == dd
->nc
[dim
] - 1) ||
492 (dir
== 1 && dd
->ci
[dim
] == 0)))
494 /* No pbc: the fist/last cell should not request atoms */
495 nsend_ptr
= nsend_zero
;
501 /* Communicate the number of indices */
502 dd_sendrecv_int(dd
,d
,dir
==0 ? dddirForward
: dddirBackward
,
503 nsend_ptr
,2,spac
->nreq
[d
][dir
],2);
504 nr
= spac
->nreq
[d
][dir
][1];
505 if (nlast
+nr
> spac
->ind_req_nalloc
)
507 spac
->ind_req_nalloc
= over_alloc_dd(nlast
+nr
);
508 srenew(spac
->ind_req
,spac
->ind_req_nalloc
);
510 /* Communicate the indices */
511 dd_sendrecv_int(dd
,d
,dir
==0 ? dddirForward
: dddirBackward
,
512 spac
->ind_req
,nsend_ptr
[1],spac
->ind_req
+nlast
,nr
);
519 fprintf(debug
,"Communicated the counts\n");
522 /* Search for the requested atoms and communicate the indices we have */
523 nat_tot_specat
= at_start
;
525 for(d
=0; d
<dd
->ndim
; d
++)
528 /* Pulse the grid forward and backward */
529 if (dd
->dim
[d
] >= dd
->npbcdim
|| dd
->nc
[dd
->dim
[d
]] > 2)
537 nat_tot_prev
= nat_tot_specat
;
538 for(dir
=ndir
-1; dir
>=0; dir
--)
540 if (nat_tot_specat
> spac
->bSendAtom_nalloc
)
542 nalloc_old
= spac
->bSendAtom_nalloc
;
543 spac
->bSendAtom_nalloc
= over_alloc_dd(nat_tot_specat
);
544 srenew(spac
->bSendAtom
,spac
->bSendAtom_nalloc
);
545 for(i
=nalloc_old
; i
<spac
->bSendAtom_nalloc
; i
++)
547 spac
->bSendAtom
[i
] = FALSE
;
550 spas
= &spac
->spas
[d
][dir
];
551 n0
= spac
->nreq
[d
][dir
][0];
552 nr
= spac
->nreq
[d
][dir
][1];
555 fprintf(debug
,"dim=%d, dir=%d, searching for %d atoms\n",
563 ireq
= spac
->ind_req
[start
+i
];
565 /* Check if this is a home atom and if so ind will be set */
566 if (!ga2la_get_home(dd
->ga2la
,ireq
,&ind
))
568 /* Search in the communicated atoms */
569 ind
= ga2la_specat
[ireq
];
573 if (i
< n0
|| !spac
->bSendAtom
[ind
])
575 if (spas
->nsend
+1 > spas
->a_nalloc
)
577 spas
->a_nalloc
= over_alloc_large(spas
->nsend
+1);
578 srenew(spas
->a
,spas
->a_nalloc
);
580 /* Store the local index so we know which coordinates
583 spas
->a
[spas
->nsend
] = ind
;
584 spac
->bSendAtom
[ind
] = TRUE
;
585 if (spas
->nsend
+1 > spac
->ibuf_nalloc
)
587 spac
->ibuf_nalloc
= over_alloc_large(spas
->nsend
+1);
588 srenew(spac
->ibuf
,spac
->ibuf_nalloc
);
590 /* Store the global index so we can send it now */
591 spac
->ibuf
[spas
->nsend
] = ireq
;
601 /* Clear the local flags */
602 for(i
=0; i
<spas
->nsend
; i
++)
604 spac
->bSendAtom
[spas
->a
[i
]] = FALSE
;
606 /* Send and receive the number of indices to communicate */
607 nsend
[1] = spas
->nsend
;
608 dd_sendrecv_int(dd
,d
,dir
==0 ? dddirBackward
: dddirForward
,
612 fprintf(debug
,"Send to node %d, %d (%d) indices, "
613 "receive from node %d, %d (%d) indices\n",
614 dd
->neighbor
[d
][1-dir
],nsend
[1],nsend
[0],
615 dd
->neighbor
[d
][dir
],buf
[1],buf
[0]);
618 for(i
=0; i
<spas
->nsend
; i
++)
620 fprintf(debug
," %d",spac
->ibuf
[i
]+1);
625 nrecv_local
+= buf
[0];
626 spas
->nrecv
= buf
[1];
627 if (nat_tot_specat
+ spas
->nrecv
> dd
->gatindex_nalloc
)
629 dd
->gatindex_nalloc
=
630 over_alloc_dd(nat_tot_specat
+ spas
->nrecv
);
631 srenew(dd
->gatindex
,dd
->gatindex_nalloc
);
633 /* Send and receive the indices */
634 dd_sendrecv_int(dd
,d
,dir
==0 ? dddirBackward
: dddirForward
,
635 spac
->ibuf
,spas
->nsend
,
636 dd
->gatindex
+nat_tot_specat
,spas
->nrecv
);
637 nat_tot_specat
+= spas
->nrecv
;
640 /* Allocate the x/f communication buffers */
641 ns
= spac
->spas
[d
][0].nsend
;
642 nr
= spac
->spas
[d
][0].nrecv
;
645 ns
+= spac
->spas
[d
][1].nsend
;
646 nr
+= spac
->spas
[d
][1].nrecv
;
648 if (vbuf_fac
*ns
> spac
->vbuf_nalloc
)
650 spac
->vbuf_nalloc
= over_alloc_dd(vbuf_fac
*ns
);
651 srenew(spac
->vbuf
,spac
->vbuf_nalloc
);
653 if (vbuf_fac
== 2 && vbuf_fac
*nr
> spac
->vbuf2_nalloc
)
655 spac
->vbuf2_nalloc
= over_alloc_dd(vbuf_fac
*nr
);
656 srenew(spac
->vbuf2
,spac
->vbuf2_nalloc
);
659 /* Make a global to local index for the communication atoms */
660 for(i
=nat_tot_prev
; i
<nat_tot_specat
; i
++)
662 ga2la_specat
[dd
->gatindex
[i
]] = i
;
666 /* Check that in the end we got the number of atoms we asked for */
667 if (nrecv_local
!= spac
->nind_req
)
671 fprintf(debug
,"Requested %d, received %d (tot recv %d)\n",
672 spac
->nind_req
,nrecv_local
,nat_tot_specat
-at_start
);
675 for(i
=0; i
<spac
->nind_req
; i
++)
677 fprintf(debug
," %s%d",
678 ga2la_specat
[spac
->ind_req
[i
]]>=0 ? "" : "!",
684 fprintf(stderr
,"\nDD cell %d %d %d: Neighboring cells do not have atoms:",
685 dd
->ci
[XX
],dd
->ci
[YY
],dd
->ci
[ZZ
]);
686 for(i
=0; i
<spac
->nind_req
; i
++)
688 if (ga2la_specat
[spac
->ind_req
[i
]] < 0)
690 fprintf(stderr
," %d",spac
->ind_req
[i
]+1);
693 fprintf(stderr
,"\n");
694 gmx_fatal(FARGS
,"DD cell %d %d %d could only obtain %d of the %d atoms that are connected via %ss from the neighboring cells. This probably means your %s lengths are too long compared to the domain decomposition cell size. Decrease the number of domain decomposition grid cells%s%s.",
695 dd
->ci
[XX
],dd
->ci
[YY
],dd
->ci
[ZZ
],
696 nrecv_local
,spac
->nind_req
,specat_type
,
698 dd
->bGridJump
? " or use the -rcon option of mdrun" : "");
701 spac
->at_start
= at_start
;
702 spac
->at_end
= nat_tot_specat
;
706 fprintf(debug
,"Done setup_specat_communication\n");
709 return nat_tot_specat
;
712 static void walk_out(int con
,int con_offset
,int a
,int offset
,int nrec
,
713 int ncon1
,const t_iatom
*ia1
,const t_iatom
*ia2
,
714 const t_blocka
*at2con
,
715 const gmx_ga2la_t ga2la
,gmx_bool bHomeConnect
,
716 gmx_domdec_constraints_t
*dc
,
717 gmx_domdec_specat_comm_t
*dcc
,
720 int a1_gl
,a2_gl
,a_loc
,i
,coni
,b
;
723 if (dc
->gc_req
[con_offset
+con
] == 0)
725 /* Add this non-home constraint to the list */
726 if (dc
->ncon
+1 > dc
->con_nalloc
)
728 dc
->con_nalloc
= over_alloc_large(dc
->ncon
+1);
729 srenew(dc
->con_gl
,dc
->con_nalloc
);
730 srenew(dc
->con_nlocat
,dc
->con_nalloc
);
732 dc
->con_gl
[dc
->ncon
] = con_offset
+ con
;
733 dc
->con_nlocat
[dc
->ncon
] = (bHomeConnect
? 1 : 0);
734 dc
->gc_req
[con_offset
+con
] = 1;
735 if (il_local
->nr
+ 3 > il_local
->nalloc
)
737 il_local
->nalloc
= over_alloc_dd(il_local
->nr
+3);
738 srenew(il_local
->iatoms
,il_local
->nalloc
);
740 iap
= constr_iatomptr(ncon1
,ia1
,ia2
,con
);
741 il_local
->iatoms
[il_local
->nr
++] = iap
[0];
742 a1_gl
= offset
+ iap
[1];
743 a2_gl
= offset
+ iap
[2];
744 /* The following indexing code can probably be optizimed */
745 if (ga2la_get_home(ga2la
,a1_gl
,&a_loc
))
747 il_local
->iatoms
[il_local
->nr
++] = a_loc
;
751 /* We set this index later */
752 il_local
->iatoms
[il_local
->nr
++] = -a1_gl
- 1;
754 if (ga2la_get_home(ga2la
,a2_gl
,&a_loc
))
756 il_local
->iatoms
[il_local
->nr
++] = a_loc
;
760 /* We set this index later */
761 il_local
->iatoms
[il_local
->nr
++] = -a2_gl
- 1;
765 /* Check to not ask for the same atom more than once */
766 if (dc
->ga2la
[offset
+a
] == -1)
769 /* Add this non-home atom to the list */
770 if (dcc
->nind_req
+1 > dcc
->ind_req_nalloc
)
772 dcc
->ind_req_nalloc
= over_alloc_large(dcc
->nind_req
+1);
773 srenew(dcc
->ind_req
,dcc
->ind_req_nalloc
);
775 dcc
->ind_req
[dcc
->nind_req
++] = offset
+ a
;
776 /* Temporarily mark with -2, we get the index later */
777 dc
->ga2la
[offset
+a
] = -2;
782 for(i
=at2con
->index
[a
]; i
<at2con
->index
[a
+1]; i
++)
788 iap
= constr_iatomptr(ncon1
,ia1
,ia2
,coni
);
797 if (!ga2la_get_home(ga2la
,offset
+b
,&a_loc
))
799 walk_out(coni
,con_offset
,b
,offset
,nrec
-1,
800 ncon1
,ia1
,ia2
,at2con
,
801 ga2la
,FALSE
,dc
,dcc
,il_local
);
808 int dd_make_local_constraints(gmx_domdec_t
*dd
,int at_start
,
810 gmx_constr_t constr
,int nrec
,
813 t_blocka
*at2con_mt
,*at2con
;
816 gmx_molblock_t
*molb
;
817 t_iatom
*ia1
,*ia2
,*iap
;
818 int nhome
,a
,a_gl
,a_mol
,a_loc
,b_lo
,offset
,mb
,molnr
,b_mol
,i
,con
,con_offset
;
819 gmx_domdec_constraints_t
*dc
;
820 int at_end
,*ga2la_specat
,j
;
822 dc
= dd
->constraints
;
824 at2con_mt
= atom2constraints_moltype(constr
);
830 if (dd
->constraint_comm
)
832 dd
->constraint_comm
->nind_req
= 0;
834 for(a
=0; a
<dd
->nat_home
; a
++)
836 a_gl
= dd
->gatindex
[a
];
838 gmx_mtop_atomnr_to_molblock_ind(mtop
,a_gl
,&mb
,&molnr
,&a_mol
);
839 molb
= &mtop
->molblock
[mb
];
841 ncon1
= mtop
->moltype
[molb
->type
].ilist
[F_CONSTR
].nr
/3;
842 ncon2
= mtop
->moltype
[molb
->type
].ilist
[F_CONSTRNC
].nr
/3;
843 if (ncon1
> 0 || ncon2
> 0)
845 ia1
= mtop
->moltype
[molb
->type
].ilist
[F_CONSTR
].iatoms
;
846 ia2
= mtop
->moltype
[molb
->type
].ilist
[F_CONSTRNC
].iatoms
;
848 /* Calculate the global constraint number offset for the molecule.
849 * This is only required for the global index to make sure
850 * that we use each constraint only once.
852 con_offset
= dc
->molb_con_offset
[mb
] + molnr
*dc
->molb_ncon_mol
[mb
];
854 /* The global atom number offset for this molecule */
855 offset
= a_gl
- a_mol
;
856 at2con
= &at2con_mt
[molb
->type
];
857 for(i
=at2con
->index
[a_mol
]; i
<at2con
->index
[a_mol
+1]; i
++)
860 iap
= constr_iatomptr(ncon1
,ia1
,ia2
,con
);
869 if (ga2la_get_home(ga2la
,offset
+b_mol
,&a_loc
))
871 /* Add this fully home constraint at the first atom */
874 if (dc
->ncon
+1 > dc
->con_nalloc
)
876 dc
->con_nalloc
= over_alloc_large(dc
->ncon
+1);
877 srenew(dc
->con_gl
,dc
->con_nalloc
);
878 srenew(dc
->con_nlocat
,dc
->con_nalloc
);
880 dc
->con_gl
[dc
->ncon
] = con_offset
+ con
;
881 dc
->con_nlocat
[dc
->ncon
] = 2;
882 if (il_local
->nr
+ 3 > il_local
->nalloc
)
884 il_local
->nalloc
= over_alloc_dd(il_local
->nr
+ 3);
885 srenew(il_local
->iatoms
,il_local
->nalloc
);
888 il_local
->iatoms
[il_local
->nr
++] = iap
[0];
889 il_local
->iatoms
[il_local
->nr
++] = (a_gl
== iap
[1] ? a
: b_lo
);
890 il_local
->iatoms
[il_local
->nr
++] = (a_gl
== iap
[1] ? b_lo
: a
);
897 /* We need the nrec constraints coupled to this constraint,
898 * so we need to walk out of the home cell by nrec+1 atoms,
899 * since already atom bg is not locally present.
900 * Therefore we call walk_out with nrec recursions to go
901 * after this first call.
903 walk_out(con
,con_offset
,b_mol
,offset
,nrec
,
904 ncon1
,ia1
,ia2
,at2con
,
905 dd
->ga2la
,TRUE
,dc
,dd
->constraint_comm
,il_local
);
914 "Constraints: home %3d border %3d atoms: %3d\n",
915 nhome
,dc
->ncon
-nhome
,
916 dd
->constraint_comm
? dd
->constraint_comm
->nind_req
: 0);
919 if (dd
->constraint_comm
) {
921 setup_specat_communication(dd
,dd
->constraint_comm
,
922 dd
->constraints
->ga2la
,
924 "constraint"," or lincs-order");
926 /* Fill in the missing indices */
927 ga2la_specat
= dd
->constraints
->ga2la
;
928 for(i
=0; i
<il_local
->nr
; i
+=3)
930 iap
= il_local
->iatoms
+ i
;
935 iap
[j
] = ga2la_specat
[-iap
[j
]-1];
948 int dd_make_local_vsites(gmx_domdec_t
*dd
,int at_start
,t_ilist
*lil
)
950 gmx_domdec_specat_comm_t
*spac
;
952 int ftype
,nral
,i
,j
,gat
,a
;
957 spac
= dd
->vsite_comm
;
958 ga2la_specat
= dd
->ga2la_vsite
;
961 /* Loop over all the home vsites */
962 for(ftype
=0; ftype
<F_NRE
; ftype
++)
964 if (interaction_function
[ftype
].flags
& IF_VSITE
)
968 for(i
=0; i
<lilf
->nr
; i
+=1+nral
)
970 iatoms
= lilf
->iatoms
+ i
;
971 /* Check if we have the other atoms */
972 for(j
=1; j
<1+nral
; j
++)
975 /* This is not a home atom,
976 * we need to ask our neighbors.
979 /* Check to not ask for the same atom more than once */
980 if (ga2la_specat
[a
] == -1)
982 /* Add this non-home atom to the list */
983 if (spac
->nind_req
+1 > spac
->ind_req_nalloc
)
985 spac
->ind_req_nalloc
=
986 over_alloc_small(spac
->nind_req
+1);
987 srenew(spac
->ind_req
,spac
->ind_req_nalloc
);
989 spac
->ind_req
[spac
->nind_req
++] = a
;
990 /* Temporarily mark with -2,
991 * we get the index later.
993 ga2la_specat
[a
] = -2;
1001 at_end
= setup_specat_communication(dd
,dd
->vsite_comm
,ga2la_specat
,
1002 at_start
,1,"vsite","");
1004 /* Fill in the missing indices */
1005 for(ftype
=0; ftype
<F_NRE
; ftype
++)
1007 if (interaction_function
[ftype
].flags
& IF_VSITE
)
1011 for(i
=0; i
<lilf
->nr
; i
+=1+nral
)
1013 iatoms
= lilf
->iatoms
+ i
;
1014 for(j
=1; j
<1+nral
; j
++)
1018 iatoms
[j
] = ga2la_specat
[-iatoms
[j
]-1];
1028 void init_domdec_constraints(gmx_domdec_t
*dd
,
1029 int natoms
,gmx_mtop_t
*mtop
,
1030 gmx_constr_t constr
)
1032 gmx_domdec_constraints_t
*dc
;
1033 gmx_molblock_t
*molb
;
1038 fprintf(debug
,"Begin init_domdec_constraints\n");
1041 snew(dd
->constraints
,1);
1042 dc
= dd
->constraints
;
1044 snew(dc
->molb_con_offset
,mtop
->nmolblock
);
1045 snew(dc
->molb_ncon_mol
,mtop
->nmolblock
);
1048 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
1050 molb
= &mtop
->molblock
[mb
];
1051 dc
->molb_con_offset
[mb
] = ncon
;
1052 dc
->molb_ncon_mol
[mb
] =
1053 mtop
->moltype
[molb
->type
].ilist
[F_CONSTR
].nr
/3 +
1054 mtop
->moltype
[molb
->type
].ilist
[F_CONSTRNC
].nr
/3;
1055 ncon
+= molb
->nmol
*dc
->molb_ncon_mol
[mb
];
1058 snew(dc
->gc_req
,ncon
);
1059 for(c
=0; c
<ncon
; c
++)
1064 snew(dc
->ga2la
,natoms
);
1065 for(a
=0; a
<natoms
; a
++)
1070 snew(dd
->constraint_comm
,1);
1073 void init_domdec_vsites(gmx_domdec_t
*dd
,int natoms
)
1076 gmx_domdec_constraints_t
*dc
;
1080 fprintf(debug
,"Begin init_domdec_vsites\n");
1083 snew(dd
->ga2la_vsite
,natoms
);
1084 for(i
=0; i
<natoms
; i
++)
1086 dd
->ga2la_vsite
[i
] = -1;
1089 snew(dd
->vsite_comm
,1);