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-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, 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.
44 #include "gromacs/domdec/ga2la.h"
45 #include "gromacs/fileio/confio.h"
46 #include "gromacs/gmxlib/network.h"
47 #include "gromacs/legacyheaders/names.h"
48 #include "gromacs/legacyheaders/types/commrec.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/mdtypes/mdatom.h"
51 #include "gromacs/pbcutil/pbc.h"
52 #include "gromacs/pulling/pull.h"
53 #include "gromacs/pulling/pull_internal.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/real.h"
57 #include "gromacs/utility/smalloc.h"
59 static void pull_reduce_real(t_commrec
*cr
,
64 if (cr
!= NULL
&& PAR(cr
))
66 if (comm
->bParticipateAll
)
68 /* Sum the contributions over all DD ranks */
74 #ifdef MPI_IN_PLACE_EXISTS
75 MPI_Allreduce(MPI_IN_PLACE
, data
, n
, GMX_MPI_REAL
, MPI_SUM
,
82 MPI_Allreduce(data
, buf
, n
, GMX_MPI_REAL
, MPI_SUM
,
85 /* Copy the result from the buffer to the input/output data */
86 for (int i
= 0; i
< n
; i
++)
93 gmx_incons("comm->bParticipateAll=FALSE without GMX_MPI");
99 static void pull_reduce_double(t_commrec
*cr
,
104 if (cr
!= NULL
&& PAR(cr
))
106 if (comm
->bParticipateAll
)
108 /* Sum the contributions over all DD ranks */
109 gmx_sumd(n
, data
, cr
);
114 #ifdef MPI_IN_PLACE_EXISTS
115 MPI_Allreduce(MPI_IN_PLACE
, data
, n
, MPI_DOUBLE
, MPI_SUM
,
122 MPI_Allreduce(data
, buf
, n
, MPI_DOUBLE
, MPI_SUM
,
125 /* Copy the result from the buffer to the input/output data */
126 for (int i
= 0; i
< n
; i
++)
133 gmx_incons("comm->bParticipateAll=FALSE without GMX_MPI");
139 static void pull_set_pbcatom(t_commrec
*cr
, pull_group_work_t
*pgrp
,
145 if (cr
!= NULL
&& DOMAINDECOMP(cr
))
147 if (ga2la_get_home(cr
->dd
->ga2la
, pgrp
->params
.pbcatom
, &a
))
149 copy_rvec(x
[a
], x_pbc
);
158 copy_rvec(x
[pgrp
->params
.pbcatom
], x_pbc
);
162 static void pull_set_pbcatoms(t_commrec
*cr
, struct pull_t
*pull
,
169 for (g
= 0; g
< pull
->ngroup
; g
++)
171 if (!pull
->group
[g
].bCalcCOM
|| pull
->group
[g
].params
.pbcatom
== -1)
173 clear_rvec(x_pbc
[g
]);
177 pull_set_pbcatom(cr
, &pull
->group
[g
], x
, x_pbc
[g
]);
182 if (cr
&& PAR(cr
) && n
> 0)
184 /* Sum over participating ranks to get x_pbc from the home ranks.
185 * This can be very expensive at high parallelization, so we only
186 * do this after each DD repartitioning.
188 pull_reduce_real(cr
, &pull
->comm
, pull
->ngroup
*DIM
, x_pbc
[0]);
192 static void make_cyl_refgrps(t_commrec
*cr
, struct pull_t
*pull
, t_mdatoms
*md
,
193 t_pbc
*pbc
, double t
, rvec
*x
)
195 /* The size and stride per coord for the reduction buffer */
196 const int stride
= 9;
197 int c
, i
, ii
, m
, start
, end
;
201 gmx_ga2la_t ga2la
= NULL
;
205 if (comm
->dbuf_cyl
== NULL
)
207 snew(comm
->dbuf_cyl
, pull
->ncoord
*stride
);
210 if (cr
&& DOMAINDECOMP(cr
))
212 ga2la
= cr
->dd
->ga2la
;
218 inv_cyl_r2
= 1/dsqr(pull
->params
.cylinder_r
);
220 /* loop over all groups to make a reference group for each*/
221 for (c
= 0; c
< pull
->ncoord
; c
++)
223 pull_coord_work_t
*pcrd
;
224 double sum_a
, wmass
, wwmass
;
225 dvec radf_fac0
, radf_fac1
;
227 pcrd
= &pull
->coord
[c
];
232 clear_dvec(radf_fac0
);
233 clear_dvec(radf_fac1
);
235 if (pcrd
->params
.eGeom
== epullgCYL
)
237 pull_group_work_t
*pref
, *pgrp
, *pdyna
;
239 /* pref will be the same group for all pull coordinates */
240 pref
= &pull
->group
[pcrd
->params
.group
[0]];
241 pgrp
= &pull
->group
[pcrd
->params
.group
[1]];
242 pdyna
= &pull
->dyna
[c
];
243 copy_rvec(pcrd
->vec
, dir
);
246 /* We calculate distances with respect to the reference location
247 * of this cylinder group (g_x), which we already have now since
248 * we reduced the other group COM over the ranks. This resolves
249 * any PBC issues and we don't need to use a PBC-atom here.
251 if (pcrd
->params
.rate
!= 0)
253 /* With rate=0, value_ref is set initially */
254 pcrd
->value_ref
= pcrd
->params
.init
+ pcrd
->params
.rate
*t
;
256 for (m
= 0; m
< DIM
; m
++)
258 g_x
[m
] = pgrp
->x
[m
] - pcrd
->vec
[m
]*pcrd
->value_ref
;
261 /* loop over all atoms in the main ref group */
262 for (i
= 0; i
< pref
->params
.nat
; i
++)
264 ii
= pref
->params
.ind
[i
];
267 if (!ga2la_get_home(ga2la
, pref
->params
.ind
[i
], &ii
))
272 if (ii
>= start
&& ii
< end
)
274 double dr2
, dr2_rel
, inp
;
277 pbc_dx_aiuc(pbc
, x
[ii
], g_x
, dx
);
278 inp
= iprod(dir
, dx
);
280 for (m
= 0; m
< DIM
; m
++)
282 /* Determine the radial components */
283 dr
[m
] = dx
[m
] - inp
*dir
[m
];
286 dr2_rel
= dr2
*inv_cyl_r2
;
290 double mass
, weight
, dweight_r
;
293 /* add to index, to sum of COM, to weight array */
294 if (pdyna
->nat_loc
>= pdyna
->nalloc_loc
)
296 pdyna
->nalloc_loc
= over_alloc_large(pdyna
->nat_loc
+1);
297 srenew(pdyna
->ind_loc
, pdyna
->nalloc_loc
);
298 srenew(pdyna
->weight_loc
, pdyna
->nalloc_loc
);
299 srenew(pdyna
->mdw
, pdyna
->nalloc_loc
);
300 srenew(pdyna
->dv
, pdyna
->nalloc_loc
);
302 pdyna
->ind_loc
[pdyna
->nat_loc
] = ii
;
304 mass
= md
->massT
[ii
];
305 /* The radial weight function is 1-2x^2+x^4,
306 * where x=r/cylinder_r. Since this function depends
307 * on the radial component, we also get radial forces
310 weight
= 1 + (-2 + dr2_rel
)*dr2_rel
;
311 dweight_r
= (-4 + 4*dr2_rel
)*inv_cyl_r2
;
312 pdyna
->weight_loc
[pdyna
->nat_loc
] = weight
;
313 sum_a
+= mass
*weight
*inp
;
314 wmass
+= mass
*weight
;
315 wwmass
+= mass
*weight
*weight
;
316 dsvmul(mass
*dweight_r
, dr
, mdw
);
317 copy_dvec(mdw
, pdyna
->mdw
[pdyna
->nat_loc
]);
318 /* Currently we only have the axial component of the
319 * distance (inp) up to an unkown offset. We add this
320 * offset after the reduction needs to determine the
321 * COM of the cylinder group.
323 pdyna
->dv
[pdyna
->nat_loc
] = inp
;
324 for (m
= 0; m
< DIM
; m
++)
326 radf_fac0
[m
] += mdw
[m
];
327 radf_fac1
[m
] += mdw
[m
]*inp
;
334 comm
->dbuf_cyl
[c
*stride
+0] = wmass
;
335 comm
->dbuf_cyl
[c
*stride
+1] = wwmass
;
336 comm
->dbuf_cyl
[c
*stride
+2] = sum_a
;
337 comm
->dbuf_cyl
[c
*stride
+3] = radf_fac0
[XX
];
338 comm
->dbuf_cyl
[c
*stride
+4] = radf_fac0
[YY
];
339 comm
->dbuf_cyl
[c
*stride
+5] = radf_fac0
[ZZ
];
340 comm
->dbuf_cyl
[c
*stride
+6] = radf_fac1
[XX
];
341 comm
->dbuf_cyl
[c
*stride
+7] = radf_fac1
[YY
];
342 comm
->dbuf_cyl
[c
*stride
+8] = radf_fac1
[ZZ
];
345 if (cr
!= NULL
&& PAR(cr
))
347 /* Sum the contributions over the ranks */
348 pull_reduce_double(cr
, comm
, pull
->ncoord
*stride
, comm
->dbuf_cyl
);
351 for (c
= 0; c
< pull
->ncoord
; c
++)
353 pull_coord_work_t
*pcrd
;
355 pcrd
= &pull
->coord
[c
];
357 if (pcrd
->params
.eGeom
== epullgCYL
)
359 pull_group_work_t
*pdyna
, *pgrp
;
360 double wmass
, wwmass
, dist
;
362 pdyna
= &pull
->dyna
[c
];
363 pgrp
= &pull
->group
[pcrd
->params
.group
[1]];
365 wmass
= comm
->dbuf_cyl
[c
*stride
+0];
366 wwmass
= comm
->dbuf_cyl
[c
*stride
+1];
367 pdyna
->mwscale
= 1.0/wmass
;
368 /* Cylinder pulling can't be used with constraints, but we set
369 * wscale and invtm anyhow, in case someone would like to use them.
371 pdyna
->wscale
= wmass
/wwmass
;
372 pdyna
->invtm
= wwmass
/(wmass
*wmass
);
374 /* We store the deviation of the COM from the reference location
375 * used above, since we need it when we apply the radial forces
376 * to the atoms in the cylinder group.
379 for (m
= 0; m
< DIM
; m
++)
381 g_x
[m
] = pgrp
->x
[m
] - pcrd
->vec
[m
]*pcrd
->value_ref
;
382 dist
= -pcrd
->vec
[m
]*comm
->dbuf_cyl
[c
*stride
+2]*pdyna
->mwscale
;
383 pdyna
->x
[m
] = g_x
[m
] - dist
;
384 pcrd
->cyl_dev
+= dist
;
386 /* Now we know the exact COM of the cylinder reference group,
387 * we can determine the radial force factor (ffrad) that when
388 * multiplied with the axial pull force will give the radial
389 * force on the pulled (non-cylinder) group.
391 for (m
= 0; m
< DIM
; m
++)
393 pcrd
->ffrad
[m
] = (comm
->dbuf_cyl
[c
*stride
+6+m
] +
394 comm
->dbuf_cyl
[c
*stride
+3+m
]*pcrd
->cyl_dev
)/wmass
;
399 fprintf(debug
, "Pull cylinder group %d:%8.3f%8.3f%8.3f m:%8.3f\n",
400 c
, pdyna
->x
[0], pdyna
->x
[1],
401 pdyna
->x
[2], 1.0/pdyna
->invtm
);
402 fprintf(debug
, "ffrad %8.3f %8.3f %8.3f\n",
403 pcrd
->ffrad
[XX
], pcrd
->ffrad
[YY
], pcrd
->ffrad
[ZZ
]);
409 static double atan2_0_2pi(double y
, double x
)
421 /* calculates center of mass of selection index from all coordinates x */
422 void pull_calc_coms(t_commrec
*cr
,
423 struct pull_t
*pull
, t_mdatoms
*md
, t_pbc
*pbc
, double t
,
432 if (comm
->rbuf
== NULL
)
434 snew(comm
->rbuf
, pull
->ngroup
);
436 if (comm
->dbuf
== NULL
)
438 snew(comm
->dbuf
, 3*pull
->ngroup
);
441 if (pull
->bRefAt
&& pull
->bSetPBCatoms
)
443 pull_set_pbcatoms(cr
, pull
, x
, comm
->rbuf
);
445 if (cr
!= NULL
&& DOMAINDECOMP(cr
))
447 /* We can keep these PBC reference coordinates fixed for nstlist
448 * steps, since atoms won't jump over PBC.
449 * This avoids a global reduction at the next nstlist-1 steps.
450 * Note that the exact values of the pbc reference coordinates
451 * are irrelevant, as long all atoms in the group are within
452 * half a box distance of the reference coordinate.
454 pull
->bSetPBCatoms
= FALSE
;
458 if (pull
->cosdim
>= 0)
462 assert(pull
->npbcdim
<= DIM
);
464 for (m
= pull
->cosdim
+1; m
< pull
->npbcdim
; m
++)
466 if (pbc
->box
[m
][pull
->cosdim
] != 0)
468 gmx_fatal(FARGS
, "Can not do cosine weighting for trilinic dimensions");
471 twopi_box
= 2.0*M_PI
/pbc
->box
[pull
->cosdim
][pull
->cosdim
];
474 for (g
= 0; g
< pull
->ngroup
; g
++)
476 pull_group_work_t
*pgrp
;
478 pgrp
= &pull
->group
[g
];
482 if (pgrp
->epgrppbc
!= epgrppbcCOS
)
485 double wmass
, wwmass
;
486 rvec x_pbc
= { 0, 0, 0 };
494 if (pgrp
->epgrppbc
== epgrppbcREFAT
)
496 /* Set the pbc atom */
497 copy_rvec(comm
->rbuf
[g
], x_pbc
);
500 for (i
= 0; i
< pgrp
->nat_loc
; i
++)
505 ii
= pgrp
->ind_loc
[i
];
506 mass
= md
->massT
[ii
];
507 if (pgrp
->weight_loc
== NULL
)
516 w
= pgrp
->weight_loc
[i
];
521 if (pgrp
->epgrppbc
== epgrppbcNONE
)
523 /* Plain COM: sum the coordinates */
524 for (m
= 0; m
< DIM
; m
++)
526 com
[m
] += wm
*x
[ii
][m
];
530 for (m
= 0; m
< DIM
; m
++)
532 comp
[m
] += wm
*xp
[ii
][m
];
540 /* Sum the difference with the reference atom */
541 pbc_dx(pbc
, x
[ii
], x_pbc
, dx
);
542 for (m
= 0; m
< DIM
; m
++)
548 /* For xp add the difference between xp and x to dx,
549 * such that we use the same periodic image,
550 * also when xp has a large displacement.
552 for (m
= 0; m
< DIM
; m
++)
554 comp
[m
] += wm
*(dx
[m
] + xp
[ii
][m
] - x
[ii
][m
]);
560 /* We do this check after the loop above to avoid more nesting.
561 * If we have a single-atom group the mass is irrelevant, so
562 * we can remove the mass factor to avoid division by zero.
563 * Note that with constraint pulling the mass does matter, but
564 * in that case a check group mass != 0 has been done before.
566 if (pgrp
->params
.nat
== 1 && pgrp
->nat_loc
== 1 && wmass
== 0)
570 /* Copy the single atom coordinate */
571 for (m
= 0; m
< DIM
; m
++)
573 com
[m
] = x
[pgrp
->ind_loc
[0]][m
];
575 /* Set all mass factors to 1 to get the correct COM */
580 if (pgrp
->weight_loc
== NULL
)
585 /* Copy local sums to a buffer for global summing */
586 copy_dvec(com
, comm
->dbuf
[g
*3]);
587 copy_dvec(comp
, comm
->dbuf
[g
*3 + 1]);
588 comm
->dbuf
[g
*3 + 2][0] = wmass
;
589 comm
->dbuf
[g
*3 + 2][1] = wwmass
;
590 comm
->dbuf
[g
*3 + 2][2] = 0;
594 /* Cosine weighting geometry */
595 double cm
, sm
, cmp
, smp
, ccm
, csm
, ssm
, csw
, snw
;
606 for (i
= 0; i
< pgrp
->nat_loc
; i
++)
611 ii
= pgrp
->ind_loc
[i
];
612 mass
= md
->massT
[ii
];
613 /* Determine cos and sin sums */
614 csw
= cos(x
[ii
][pull
->cosdim
]*twopi_box
);
615 snw
= sin(x
[ii
][pull
->cosdim
]*twopi_box
);
624 csw
= cos(xp
[ii
][pull
->cosdim
]*twopi_box
);
625 snw
= sin(xp
[ii
][pull
->cosdim
]*twopi_box
);
631 /* Copy local sums to a buffer for global summing */
632 comm
->dbuf
[g
*3 ][0] = cm
;
633 comm
->dbuf
[g
*3 ][1] = sm
;
634 comm
->dbuf
[g
*3 ][2] = 0;
635 comm
->dbuf
[g
*3+1][0] = ccm
;
636 comm
->dbuf
[g
*3+1][1] = csm
;
637 comm
->dbuf
[g
*3+1][2] = ssm
;
638 comm
->dbuf
[g
*3+2][0] = cmp
;
639 comm
->dbuf
[g
*3+2][1] = smp
;
640 comm
->dbuf
[g
*3+2][2] = 0;
645 pull_reduce_double(cr
, comm
, pull
->ngroup
*3*DIM
, comm
->dbuf
[0]);
647 for (g
= 0; g
< pull
->ngroup
; g
++)
649 pull_group_work_t
*pgrp
;
651 pgrp
= &pull
->group
[g
];
652 if (pgrp
->params
.nat
> 0 && pgrp
->bCalcCOM
)
654 if (pgrp
->epgrppbc
!= epgrppbcCOS
)
656 double wmass
, wwmass
;
659 /* Determine the inverse mass */
660 wmass
= comm
->dbuf
[g
*3+2][0];
661 wwmass
= comm
->dbuf
[g
*3+2][1];
662 pgrp
->mwscale
= 1.0/wmass
;
663 /* invtm==0 signals a frozen group, so then we should keep it zero */
664 if (pgrp
->invtm
!= 0)
666 pgrp
->wscale
= wmass
/wwmass
;
667 pgrp
->invtm
= wwmass
/(wmass
*wmass
);
669 /* Divide by the total mass */
670 for (m
= 0; m
< DIM
; m
++)
672 pgrp
->x
[m
] = comm
->dbuf
[g
*3 ][m
]*pgrp
->mwscale
;
675 pgrp
->xp
[m
] = comm
->dbuf
[g
*3+1][m
]*pgrp
->mwscale
;
677 if (pgrp
->epgrppbc
== epgrppbcREFAT
)
679 pgrp
->x
[m
] += comm
->rbuf
[g
][m
];
682 pgrp
->xp
[m
] += comm
->rbuf
[g
][m
];
689 /* Cosine weighting geometry */
690 double csw
, snw
, wmass
, wwmass
;
693 /* Determine the optimal location of the cosine weight */
694 csw
= comm
->dbuf
[g
*3][0];
695 snw
= comm
->dbuf
[g
*3][1];
696 pgrp
->x
[pull
->cosdim
] = atan2_0_2pi(snw
, csw
)/twopi_box
;
697 /* Set the weights for the local atoms */
698 wmass
= sqrt(csw
*csw
+ snw
*snw
);
699 wwmass
= (comm
->dbuf
[g
*3+1][0]*csw
*csw
+
700 comm
->dbuf
[g
*3+1][1]*csw
*snw
+
701 comm
->dbuf
[g
*3+1][2]*snw
*snw
)/(wmass
*wmass
);
703 pgrp
->mwscale
= 1.0/wmass
;
704 pgrp
->wscale
= wmass
/wwmass
;
705 pgrp
->invtm
= wwmass
/(wmass
*wmass
);
706 /* Set the weights for the local atoms */
709 for (i
= 0; i
< pgrp
->nat_loc
; i
++)
711 ii
= pgrp
->ind_loc
[i
];
712 pgrp
->weight_loc
[i
] = csw
*cos(twopi_box
*x
[ii
][pull
->cosdim
]) +
713 snw
*sin(twopi_box
*x
[ii
][pull
->cosdim
]);
717 csw
= comm
->dbuf
[g
*3+2][0];
718 snw
= comm
->dbuf
[g
*3+2][1];
719 pgrp
->xp
[pull
->cosdim
] = atan2_0_2pi(snw
, csw
)/twopi_box
;
724 fprintf(debug
, "Pull group %d wmass %f invtm %f\n",
725 g
, 1.0/pgrp
->mwscale
, pgrp
->invtm
);
732 /* Calculate the COMs for the cyclinder reference groups */
733 make_cyl_refgrps(cr
, pull
, md
, pbc
, t
, x
);