Moved mdatom.h from legacyheader/types to mdtypes.
[gromacs.git] / src / gromacs / pulling / pullutil.cpp
blob4ea120122057ff949c122558269ac608777117a4
1 /*
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.
37 #include "gmxpre.h"
39 #include "config.h"
41 #include <assert.h>
42 #include <stdlib.h>
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,
60 pull_comm_t *comm,
61 int n,
62 real *data)
64 if (cr != NULL && PAR(cr))
66 if (comm->bParticipateAll)
68 /* Sum the contributions over all DD ranks */
69 gmx_sum(n, data, cr);
71 else
73 #ifdef GMX_MPI
74 #ifdef MPI_IN_PLACE_EXISTS
75 MPI_Allreduce(MPI_IN_PLACE, data, n, GMX_MPI_REAL, MPI_SUM,
76 comm->mpi_comm_com);
77 #else
78 real *buf;
80 snew(buf, n);
82 MPI_Allreduce(data, buf, n, GMX_MPI_REAL, MPI_SUM,
83 comm->mpi_comm_com);
85 /* Copy the result from the buffer to the input/output data */
86 for (int i = 0; i < n; i++)
88 data[i] = buf[i];
90 sfree(buf);
91 #endif
92 #else
93 gmx_incons("comm->bParticipateAll=FALSE without GMX_MPI");
94 #endif
99 static void pull_reduce_double(t_commrec *cr,
100 pull_comm_t *comm,
101 int n,
102 double *data)
104 if (cr != NULL && PAR(cr))
106 if (comm->bParticipateAll)
108 /* Sum the contributions over all DD ranks */
109 gmx_sumd(n, data, cr);
111 else
113 #ifdef GMX_MPI
114 #ifdef MPI_IN_PLACE_EXISTS
115 MPI_Allreduce(MPI_IN_PLACE, data, n, MPI_DOUBLE, MPI_SUM,
116 comm->mpi_comm_com);
117 #else
118 double *buf;
120 snew(buf, n);
122 MPI_Allreduce(data, buf, n, MPI_DOUBLE, MPI_SUM,
123 comm->mpi_comm_com);
125 /* Copy the result from the buffer to the input/output data */
126 for (int i = 0; i < n; i++)
128 data[i] = buf[i];
130 sfree(buf);
131 #endif
132 #else
133 gmx_incons("comm->bParticipateAll=FALSE without GMX_MPI");
134 #endif
139 static void pull_set_pbcatom(t_commrec *cr, pull_group_work_t *pgrp,
140 rvec *x,
141 rvec x_pbc)
143 int a;
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);
151 else
153 clear_rvec(x_pbc);
156 else
158 copy_rvec(x[pgrp->params.pbcatom], x_pbc);
162 static void pull_set_pbcatoms(t_commrec *cr, struct pull_t *pull,
163 rvec *x,
164 rvec *x_pbc)
166 int g, n;
168 n = 0;
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]);
175 else
177 pull_set_pbcatom(cr, &pull->group[g], x, x_pbc[g]);
178 n++;
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;
198 rvec g_x, dx, dir;
199 double inv_cyl_r2;
200 pull_comm_t *comm;
201 gmx_ga2la_t ga2la = NULL;
203 comm = &pull->comm;
205 if (comm->dbuf_cyl == NULL)
207 snew(comm->dbuf_cyl, pull->ncoord*stride);
210 if (cr && DOMAINDECOMP(cr))
212 ga2la = cr->dd->ga2la;
215 start = 0;
216 end = md->homenr;
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];
229 sum_a = 0;
230 wmass = 0;
231 wwmass = 0;
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);
244 pdyna->nat_loc = 0;
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];
265 if (ga2la)
267 if (!ga2la_get_home(ga2la, pref->params.ind[i], &ii))
269 ii = -1;
272 if (ii >= start && ii < end)
274 double dr2, dr2_rel, inp;
275 dvec dr;
277 pbc_dx_aiuc(pbc, x[ii], g_x, dx);
278 inp = iprod(dir, dx);
279 dr2 = 0;
280 for (m = 0; m < DIM; m++)
282 /* Determine the radial components */
283 dr[m] = dx[m] - inp*dir[m];
284 dr2 += dr[m]*dr[m];
286 dr2_rel = dr2*inv_cyl_r2;
288 if (dr2_rel < 1)
290 double mass, weight, dweight_r;
291 dvec mdw;
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
308 * on both groups.
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;
329 pdyna->nat_loc++;
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.
378 pcrd->cyl_dev = 0;
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;
397 if (debug)
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)
411 double a;
413 a = atan2(y, x);
414 if (a < 0)
416 a += 2.0*M_PI;
418 return a;
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,
424 rvec x[], rvec *xp)
426 int g;
427 real twopi_box = 0;
428 pull_comm_t *comm;
430 comm = &pull->comm;
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)
460 int m;
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];
480 if (pgrp->bCalcCOM)
482 if (pgrp->epgrppbc != epgrppbcCOS)
484 dvec com, comp;
485 double wmass, wwmass;
486 rvec x_pbc = { 0, 0, 0 };
487 int i;
489 clear_dvec(com);
490 clear_dvec(comp);
491 wmass = 0;
492 wwmass = 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++)
502 int ii, m;
503 real mass, wm;
505 ii = pgrp->ind_loc[i];
506 mass = md->massT[ii];
507 if (pgrp->weight_loc == NULL)
509 wm = mass;
510 wmass += wm;
512 else
514 real w;
516 w = pgrp->weight_loc[i];
517 wm = w*mass;
518 wmass += wm;
519 wwmass += wm*w;
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];
528 if (xp)
530 for (m = 0; m < DIM; m++)
532 comp[m] += wm*xp[ii][m];
536 else
538 rvec dx;
540 /* Sum the difference with the reference atom */
541 pbc_dx(pbc, x[ii], x_pbc, dx);
542 for (m = 0; m < DIM; m++)
544 com[m] += wm*dx[m];
546 if (xp)
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)
568 int m;
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 */
576 wmass = 1;
577 wwmass = 1;
580 if (pgrp->weight_loc == NULL)
582 wwmass = wmass;
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;
592 else
594 /* Cosine weighting geometry */
595 double cm, sm, cmp, smp, ccm, csm, ssm, csw, snw;
596 int i;
598 cm = 0;
599 sm = 0;
600 cmp = 0;
601 smp = 0;
602 ccm = 0;
603 csm = 0;
604 ssm = 0;
606 for (i = 0; i < pgrp->nat_loc; i++)
608 int ii;
609 real mass;
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);
616 cm += csw*mass;
617 sm += snw*mass;
618 ccm += csw*csw*mass;
619 csm += csw*snw*mass;
620 ssm += snw*snw*mass;
622 if (xp)
624 csw = cos(xp[ii][pull->cosdim]*twopi_box);
625 snw = sin(xp[ii][pull->cosdim]*twopi_box);
626 cmp += csw*mass;
627 smp += snw*mass;
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;
657 int m;
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;
673 if (xp)
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];
680 if (xp)
682 pgrp->xp[m] += comm->rbuf[g][m];
687 else
689 /* Cosine weighting geometry */
690 double csw, snw, wmass, wwmass;
691 int i, ii;
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 */
707 csw *= pgrp->invtm;
708 snw *= pgrp->invtm;
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]);
715 if (xp)
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;
722 if (debug)
724 fprintf(debug, "Pull group %d wmass %f invtm %f\n",
725 g, 1.0/pgrp->mwscale, pgrp->invtm);
730 if (pull->bCylinder)
732 /* Calculate the COMs for the cyclinder reference groups */
733 make_cyl_refgrps(cr, pull, md, pbc, t, x);