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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
52 #include "gmx_fatal.h"
60 #include "gmx_ga2la.h"
62 static void pull_set_pbcatom(t_commrec *cr, t_pullgrp *pg,
63 t_mdatoms *md, rvec *x,
69 if (DOMAINDECOMP(cr)) {
70 if (!ga2la_get_home(cr->dd->ga2la,pg->pbcatom,&a)) {
77 if (a >= md->start && a < md->start+md->homenr) {
78 copy_rvec(x[a],x_pbc);
83 copy_rvec(x[pg->pbcatom],x_pbc);
87 static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull,
88 t_mdatoms *md, rvec *x,
94 for(g=0; g<1+pull->ngrp; g++) {
95 if ((g==0 && PULL_CYL(pull)) || pull->grp[g].pbcatom == -1) {
98 pull_set_pbcatom(cr,&pull->grp[g],md,x,x_pbc[g]);
99 for(m=0; m<DIM; m++) {
100 if (pull->dim[m] == 0) {
108 if (cr && PAR(cr) && n > 0) {
109 /* Sum over the nodes to get x_pbc from the home node of pbcatom */
110 gmx_sum((1+pull->ngrp)*DIM,x_pbc[0],cr);
114 /* switch function, x between r and w */
115 static real get_weight(real x, real r1, real r0)
124 weight = (r0 - x)/(r0 - r1);
129 static void make_cyl_refgrps(t_commrec *cr,t_pull *pull,t_mdatoms *md,
130 t_pbc *pbc,double t,rvec *x,rvec *xp)
132 int g,i,ii,m,start,end;
134 double r0_2,sum_a,sum_ap,dr2,mass,weight,wmass,wwmass,inp;
135 t_pullgrp *pref,*pgrp,*pdyna;
136 gmx_ga2la_t ga2la=NULL;
138 if (pull->dbuf_cyl == NULL) {
139 snew(pull->dbuf_cyl,pull->ngrp*4);
142 if (cr && DOMAINDECOMP(cr))
143 ga2la = cr->dd->ga2la;
146 end = md->homenr + start;
148 r0_2 = dsqr(pull->cyl_r0);
150 /* loop over all groups to make a reference group for each*/
151 pref = &pull->grp[0];
152 for(g=1; g<1+pull->ngrp; g++) {
153 pgrp = &pull->grp[g];
154 pdyna = &pull->dyna[g];
155 copy_rvec(pgrp->vec,dir);
162 for(m=0; m<DIM; m++) {
163 g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t);
166 /* loop over all atoms in the main ref group */
167 for(i=0; i<pref->nat; i++) {
168 ii = pull->grp[0].ind[i];
170 if (!ga2la_get_home(ga2la,pref->ind[i],&ii)) {
174 if (ii >= start && ii < end) {
175 pbc_dx_aiuc(pbc,x[ii],g_x,dx);
178 for(m=0; m<DIM; m++) {
179 dr2 += dsqr(dx[m] - inp*dir[m]);
183 /* add to index, to sum of COM, to weight array */
184 if (pdyna->nat_loc >= pdyna->nalloc_loc) {
185 pdyna->nalloc_loc = over_alloc_large(pdyna->nat_loc+1);
186 srenew(pdyna->ind_loc,pdyna->nalloc_loc);
187 srenew(pdyna->weight_loc,pdyna->nalloc_loc);
189 pdyna->ind_loc[pdyna->nat_loc] = ii;
190 mass = md->massT[ii];
191 weight = get_weight(sqrt(dr2),pull->cyl_r1,pull->cyl_r0);
192 pdyna->weight_loc[pdyna->nat_loc] = weight;
193 sum_a += mass*weight*inp;
195 pbc_dx_aiuc(pbc,xp[ii],g_x,dx);
197 sum_ap += mass*weight*inp;
199 wmass += mass*weight;
200 wwmass += mass*sqr(weight);
205 pull->dbuf_cyl[(g-1)*4+0] = wmass;
206 pull->dbuf_cyl[(g-1)*4+1] = wwmass;
207 pull->dbuf_cyl[(g-1)*4+2] = sum_a;
208 pull->dbuf_cyl[(g-1)*4+3] = sum_ap;
212 /* Sum the contributions over the nodes */
213 gmx_sumd(pull->ngrp*4,pull->dbuf_cyl,cr);
216 for(g=1; g<1+pull->ngrp; g++) {
217 pgrp = &pull->grp[g];
218 pdyna = &pull->dyna[g];
220 wmass = pull->dbuf_cyl[(g-1)*4+0];
221 wwmass = pull->dbuf_cyl[(g-1)*4+1];
222 pdyna->wscale = wmass/wwmass;
223 pdyna->invtm = 1.0/(pdyna->wscale*wmass);
225 for(m=0; m<DIM; m++) {
226 g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t);
227 pdyna->x[m] = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+2]/wmass;
229 pdyna->xp[m] = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+3]/wmass;
234 fprintf(debug,"Pull cylinder group %d:%8.3f%8.3f%8.3f m:%8.3f\n",
235 g,pdyna->x[0],pdyna->x[1],
236 pdyna->x[2],1.0/pdyna->invtm);
241 static double atan2_0_2pi(double y,double x)
252 /* calculates center of mass of selection index from all coordinates x */
253 void pull_calc_coms(t_commrec *cr,
254 t_pull *pull, t_mdatoms *md, t_pbc *pbc,double t,
258 real mass,w,wm,twopi_box=0;
259 double wmass,wwmass,invwmass;
261 double cm,sm,cmp,smp,ccm,csm,ssm,csw,snw;
262 rvec *xx[2],x_pbc={0,0,0},dx;
265 if (pull->rbuf == NULL) {
266 snew(pull->rbuf,1+pull->ngrp);
268 if (pull->dbuf == NULL) {
269 snew(pull->dbuf,3*(1+pull->ngrp));
273 pull_set_pbcatoms(cr,pull,md,x,pull->rbuf);
276 if (pull->cosdim >= 0) {
277 for(m=pull->cosdim+1; m<pull->npbcdim; m++) {
278 if (pbc->box[m][pull->cosdim] != 0) {
279 gmx_fatal(FARGS,"Can not do cosine weighting for trilinic dimensions");
282 twopi_box = 2.0*M_PI/pbc->box[pull->cosdim][pull->cosdim];
285 for (g=0; g<1+pull->ngrp; g++) {
286 pgrp = &pull->grp[g];
298 if (!(g==0 && PULL_CYL(pull))) {
299 if (pgrp->epgrppbc == epgrppbcREFAT) {
300 /* Set the pbc atom */
301 copy_rvec(pull->rbuf[g],x_pbc);
304 for(i=0; i<pgrp->nat_loc; i++) {
305 ii = pgrp->ind_loc[i];
306 mass = md->massT[ii];
307 if (pgrp->epgrppbc != epgrppbcCOS) {
308 if (pgrp->weight_loc) {
309 w = pgrp->weight_loc[i];
314 if (pgrp->epgrppbc == epgrppbcNONE) {
315 /* Plain COM: sum the coordinates */
317 com[m] += wm*x[ii][m];
320 comp[m] += wm*xp[ii][m];
323 /* Sum the difference with the reference atom */
324 pbc_dx(pbc,x[ii],x_pbc,dx);
325 for(m=0; m<DIM; m++) {
329 /* For xp add the difference between xp and x to dx,
330 * such that we use the same periodic image,
331 * also when xp has a large displacement.
333 for(m=0; m<DIM; m++) {
334 comp[m] += wm*(dx[m] + xp[ii][m] - x[ii][m]);
339 /* Determine cos and sin sums */
340 csw = cos(x[ii][pull->cosdim]*twopi_box);
341 snw = sin(x[ii][pull->cosdim]*twopi_box);
349 csw = cos(xp[ii][pull->cosdim]*twopi_box);
350 snw = sin(xp[ii][pull->cosdim]*twopi_box);
358 /* Copy local sums to a buffer for global summing */
359 switch (pgrp->epgrppbc) {
362 copy_dvec(com,pull->dbuf[g*3]);
363 copy_dvec(comp,pull->dbuf[g*3+1]);
364 pull->dbuf[g*3+2][0] = wmass;
365 pull->dbuf[g*3+2][1] = wwmass;
366 pull->dbuf[g*3+2][2] = 0;
369 pull->dbuf[g*3 ][0] = cm;
370 pull->dbuf[g*3 ][1] = sm;
371 pull->dbuf[g*3 ][2] = 0;
372 pull->dbuf[g*3+1][0] = ccm;
373 pull->dbuf[g*3+1][1] = csm;
374 pull->dbuf[g*3+1][2] = ssm;
375 pull->dbuf[g*3+2][0] = cmp;
376 pull->dbuf[g*3+2][1] = smp;
377 pull->dbuf[g*3+2][2] = 0;
383 /* Sum the contributions over the nodes */
384 gmx_sumd((1+pull->ngrp)*3*DIM,pull->dbuf[0],cr);
387 for (g=0; g<1+pull->ngrp; g++) {
388 pgrp = &pull->grp[g];
389 if (pgrp->nat > 0 && !(g==0 && PULL_CYL(pull))) {
390 if (pgrp->epgrppbc != epgrppbcCOS) {
391 /* Determine the inverse mass */
392 wmass = pull->dbuf[g*3+2][0];
393 wwmass = pull->dbuf[g*3+2][1];
395 /* invtm==0 signals a frozen group, so then we should keep it zero */
396 if (pgrp->invtm > 0) {
397 pgrp->wscale = wmass/wwmass;
398 pgrp->invtm = 1.0/(pgrp->wscale*wmass);
400 /* Divide by the total mass */
401 for(m=0; m<DIM; m++) {
402 pgrp->x[m] = pull->dbuf[g*3 ][m]*invwmass;
404 pgrp->xp[m] = pull->dbuf[g*3+1][m]*invwmass;
406 if (pgrp->epgrppbc == epgrppbcREFAT) {
407 pgrp->x[m] += pull->rbuf[g][m];
409 pgrp->xp[m] += pull->rbuf[g][m];
414 /* Determine the optimal location of the cosine weight */
415 csw = pull->dbuf[g*3][0];
416 snw = pull->dbuf[g*3][1];
417 pgrp->x[pull->cosdim] = atan2_0_2pi(snw,csw)/twopi_box;
418 /* Set the weights for the local atoms */
419 wmass = sqrt(csw*csw + snw*snw);
420 wwmass = (pull->dbuf[g*3+1][0]*csw*csw +
421 pull->dbuf[g*3+1][1]*csw*snw +
422 pull->dbuf[g*3+1][2]*snw*snw)/(wmass*wmass);
423 pgrp->wscale = wmass/wwmass;
424 pgrp->invtm = 1.0/(pgrp->wscale*wmass);
425 /* Set the weights for the local atoms */
428 for(i=0; i<pgrp->nat_loc; i++) {
429 ii = pgrp->ind_loc[i];
430 pgrp->weight_loc[i] = csw*cos(twopi_box*x[ii][pull->cosdim]) +
431 snw*sin(twopi_box*x[ii][pull->cosdim]);
434 csw = pull->dbuf[g*3+2][0];
435 snw = pull->dbuf[g*3+2][1];
436 pgrp->xp[pull->cosdim] = atan2_0_2pi(snw,csw)/twopi_box;
440 fprintf(debug,"Pull group %d wmass %f wwmass %f invtm %f\n",
441 g,wmass,wwmass,pgrp->invtm);
446 if (PULL_CYL(pull)) {
447 /* Calculate the COMs for the cyclinder reference groups */
448 make_cyl_refgrps(cr,pull,md,pbc,t,x,xp);