2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2012,2014,2015,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "groupcoord.h"
41 #include "gromacs/domdec/ga2la.h"
42 #include "gromacs/gmxlib/network.h"
43 #include "gromacs/math/vec.h"
44 #include "gromacs/mdtypes/commrec.h"
45 #include "gromacs/pbcutil/pbc.h"
46 #include "gromacs/utility/smalloc.h"
48 #define MIN(a, b) (((a) < (b)) ? (a) : (b))
52 /* Select the indices of the group's atoms which are local and store them in
53 * anrs_loc[0..nr_loc]. The indices are saved in coll_ind[] for later reduction
54 * in communicate_group_positions()
56 void dd_make_local_group_indices(gmx_ga2la_t
*ga2la
,
57 const int nr
, /* IN: Total number of atoms in the group */
58 int anrs
[], /* IN: Global atom numbers of the groups atoms */
59 int *nr_loc
, /* OUT: Number of group atoms found locally */
60 int *anrs_loc
[], /* OUT: Local atom numbers of the group */
61 int *nalloc_loc
, /* IN+OUT: Allocation size of anrs_loc */
62 int coll_ind
[]) /* OUT (opt): Where is this position found in the collective array? */
68 /* Loop over all the atom indices of the group to check
69 * which ones are on the local node */
71 for (i
= 0; i
< nr
; i
++)
73 if (ga2la_get_home(ga2la
, anrs
[i
], &ii
))
75 /* The atom with this index is a home atom */
76 if (localnr
>= *nalloc_loc
) /* Check whether memory suffices */
78 *nalloc_loc
= over_alloc_dd(localnr
+1);
79 /* We never need more memory than the number of atoms in the group */
80 *nalloc_loc
= MIN(*nalloc_loc
, nr
);
81 srenew(*anrs_loc
, *nalloc_loc
);
83 /* Save the atoms index in the local atom numbers array */
84 (*anrs_loc
)[localnr
] = ii
;
86 if (coll_ind
!= nullptr)
88 /* Keep track of where this local atom belongs in the collective index array.
89 * This is needed when reducing the local arrays to a collective/global array
90 * in communicate_group_positions */
91 coll_ind
[localnr
] = i
;
94 /* add one to the local atom count */
99 /* Return the number of local atoms that were found */
104 static void get_shifts_group(
107 rvec
*xcoll
, /* IN: Collective set of positions [0..nr] */
108 int nr
, /* IN: Total number of atoms in the group */
109 rvec
*xcoll_old
, /* IN: Positions from the last time step [0...nr] */
110 ivec
*shifts
) /* OUT: Shifts for xcoll */
116 /* Get the shifts such that each atom is within closest
117 * distance to its position at the last NS time step after shifting.
118 * If we start with a whole group, and always keep track of
119 * shift changes, the group will stay whole this way */
120 for (i
= 0; i
< nr
; i
++)
122 clear_ivec(shifts
[i
]);
125 for (i
= 0; i
< nr
; i
++)
127 /* The distance this atom moved since the last time step */
128 /* If this is more than just a bit, it has changed its home pbc box */
129 rvec_sub(xcoll
[i
], xcoll_old
[i
], dx
);
131 for (m
= npbcdim
-1; m
>= 0; m
--)
133 while (dx
[m
] < -0.5*box
[m
][m
])
135 for (d
= 0; d
< DIM
; d
++)
141 while (dx
[m
] >= 0.5*box
[m
][m
])
143 for (d
= 0; d
< DIM
; d
++)
154 static void shift_positions_group(
156 rvec x
[], /* The positions [0..nr] */
157 ivec
*is
, /* The shifts [0..nr] */
158 int nr
) /* The number of positions and shifts */
163 /* Loop over the group's atoms */
166 for (i
= 0; i
< nr
; i
++)
172 x
[i
][XX
] = x
[i
][XX
]+tx
*box
[XX
][XX
]+ty
*box
[YY
][XX
]+tz
*box
[ZZ
][XX
];
173 x
[i
][YY
] = x
[i
][YY
]+ty
*box
[YY
][YY
]+tz
*box
[ZZ
][YY
];
174 x
[i
][ZZ
] = x
[i
][ZZ
]+tz
*box
[ZZ
][ZZ
];
179 for (i
= 0; i
< nr
; i
++)
185 x
[i
][XX
] = x
[i
][XX
]+tx
*box
[XX
][XX
];
186 x
[i
][YY
] = x
[i
][YY
]+ty
*box
[YY
][YY
];
187 x
[i
][ZZ
] = x
[i
][ZZ
]+tz
*box
[ZZ
][ZZ
];
193 /* Assemble the positions of the group such that every node has all of them.
194 * The atom indices are retrieved from anrs_loc[0..nr_loc]
195 * Note that coll_ind[i] = i is needed in the serial case */
196 extern void communicate_group_positions(
197 const t_commrec
*cr
, /* Pointer to MPI communication data */
198 rvec
*xcoll
, /* Collective array of positions */
199 ivec
*shifts
, /* Collective array of shifts for xcoll (can be NULL) */
200 ivec
*extra_shifts
, /* (optional) Extra shifts since last time step */
201 const gmx_bool bNS
, /* (optional) NS step, the shifts have changed */
202 const rvec
*x_loc
, /* Local positions on this node */
203 const int nr
, /* Total number of atoms in the group */
204 const int nr_loc
, /* Local number of atoms in the group */
205 int *anrs_loc
, /* Local atom numbers */
206 int *coll_ind
, /* Collective index */
207 rvec
*xcoll_old
, /* (optional) Positions from the last time step,
208 used to make group whole */
209 matrix box
) /* (optional) The box */
214 /* Zero out the groups' global position array */
215 clear_rvecs(nr
, xcoll
);
217 /* Put the local positions that this node has into the right place of
218 * the collective array. Note that in the serial case, coll_ind[i] = i */
219 for (i
= 0; i
< nr_loc
; i
++)
221 copy_rvec(x_loc
[anrs_loc
[i
]], xcoll
[coll_ind
[i
]]);
226 /* Add the arrays from all nodes together */
227 gmx_sum(nr
*3, xcoll
[0], cr
);
229 /* Now we have all the positions of the group in the xcoll array present on all
232 * The rest of the code is for making the group whole again in case atoms changed
233 * their PBC representation / crossed a box boundary. We only do that if the
234 * shifts array is allocated. */
235 if (nullptr != shifts
)
237 /* To make the group whole, start with a whole group and each
238 * step move the assembled positions at closest distance to the positions
239 * from the last step. First shift the positions with the saved shift
240 * vectors (these are 0 when this routine is called for the first time!) */
241 shift_positions_group(box
, xcoll
, shifts
, nr
);
243 /* Now check if some shifts changed since the last step.
244 * This only needs to be done when the shifts are expected to have changed,
245 * i.e. after neighbor searching */
248 get_shifts_group(3, box
, xcoll
, nr
, xcoll_old
, extra_shifts
);
250 /* Shift with the additional shifts such that we get a whole group now */
251 shift_positions_group(box
, xcoll
, extra_shifts
, nr
);
253 /* Add the shift vectors together for the next time step */
254 for (i
= 0; i
< nr
; i
++)
256 shifts
[i
][XX
] += extra_shifts
[i
][XX
];
257 shifts
[i
][YY
] += extra_shifts
[i
][YY
];
258 shifts
[i
][ZZ
] += extra_shifts
[i
][ZZ
];
261 /* Store current correctly-shifted positions for comparison in the next NS time step */
262 for (i
= 0; i
< nr
; i
++)
264 copy_rvec(xcoll
[i
], xcoll_old
[i
]);
271 /* Determine the (weighted) sum vector from positions x */
272 extern double get_sum_of_positions(rvec x
[], real weight
[], const int nat
, dvec dsumvec
)
276 double weight_sum
= 0.0;
279 /* Zero out the center */
282 /* Loop over all atoms and add their weighted position vectors */
283 if (weight
!= nullptr)
285 for (i
= 0; i
< nat
; i
++)
287 weight_sum
+= weight
[i
];
288 svmul(weight
[i
], x
[i
], x_weighted
);
289 dsumvec
[XX
] += x_weighted
[XX
];
290 dsumvec
[YY
] += x_weighted
[YY
];
291 dsumvec
[ZZ
] += x_weighted
[ZZ
];
296 for (i
= 0; i
< nat
; i
++)
298 dsumvec
[XX
] += x
[i
][XX
];
299 dsumvec
[YY
] += x
[i
][YY
];
300 dsumvec
[ZZ
] += x
[i
][ZZ
];
307 /* Determine center of structure from collective positions x */
308 extern void get_center(rvec x
[], real weight
[], const int nr
, rvec rcenter
)
311 double weight_sum
, denom
;
314 weight_sum
= get_sum_of_positions(x
, weight
, nr
, dcenter
);
316 if (weight
!= nullptr)
318 denom
= weight_sum
; /* Divide by the sum of weight */
322 denom
= nr
; /* Divide by the number of atoms */
325 dsvmul(1.0/denom
, dcenter
, dcenter
);
327 rcenter
[XX
] = dcenter
[XX
];
328 rcenter
[YY
] = dcenter
[YY
];
329 rcenter
[ZZ
] = dcenter
[ZZ
];
333 /* Get the center from local positions that already have the correct
334 * PBC representation */
335 extern void get_center_comm(
337 rvec x_loc
[], /* Local positions */
338 real weight_loc
[], /* Local masses or other weights */
339 int nr_loc
, /* Local number of atoms */
340 int nr_group
, /* Total number of atoms of the group */
341 rvec center
) /* Weighted center */
343 double weight_sum
, denom
;
348 weight_sum
= get_sum_of_positions(x_loc
, weight_loc
, nr_loc
, dsumvec
);
350 /* Add the local contributions from all nodes. Put the sum vector and the
351 * weight in a buffer array so that we get along with a single communication
355 buf
[0] = dsumvec
[XX
];
356 buf
[1] = dsumvec
[YY
];
357 buf
[2] = dsumvec
[ZZ
];
360 /* Communicate buffer */
361 gmx_sumd(4, buf
, cr
);
363 dsumvec
[XX
] = buf
[0];
364 dsumvec
[YY
] = buf
[1];
365 dsumvec
[ZZ
] = buf
[2];
369 if (weight_loc
!= nullptr)
371 denom
= 1.0/weight_sum
; /* Divide by the sum of weight to get center of mass e.g. */
375 denom
= 1.0/nr_group
; /* Divide by the number of atoms to get the geometrical center */
378 center
[XX
] = dsumvec
[XX
]*denom
;
379 center
[YY
] = dsumvec
[YY
]*denom
;
380 center
[ZZ
] = dsumvec
[ZZ
]*denom
;
384 /* Translate x with transvec */
385 extern void translate_x(rvec x
[], const int nr
, const rvec transvec
)
390 for (i
= 0; i
< nr
; i
++)
392 rvec_inc(x
[i
], transvec
);
397 extern void rotate_x(rvec x
[], const int nr
, matrix rmat
)
403 /* Apply the rotation matrix */
404 for (i
= 0; i
< nr
; i
++)
406 for (j
= 0; j
< 3; j
++)
410 for (j
= 0; j
< 3; j
++)
413 for (k
= 0; k
< 3; k
++)
415 x
[i
][j
] += rmat
[k
][j
]*x_old
[k
];