1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
59 #include "mtop_util.h"
61 #include "gmx_ga2la.h"
65 static void pull_print_x_grp(FILE *out,gmx_bool bRef,ivec dim,t_pullgrp *pgrp)
73 fprintf(out,"\t%g",bRef ? pgrp->x[m] : pgrp->dr[m]);
78 static void pull_print_x(FILE *out,t_pull *pull,double t)
82 fprintf(out, "%.4f", t);
86 for (g=1; g<1+pull->ngrp; g++)
88 pull_print_x_grp(out,TRUE ,pull->dim,&pull->dyna[g]);
89 pull_print_x_grp(out,FALSE,pull->dim,&pull->grp[g]);
94 for (g=0; g<1+pull->ngrp; g++)
96 if (pull->grp[g].nat > 0)
98 pull_print_x_grp(out,g==0,pull->dim,&pull->grp[g]);
105 static void pull_print_f(FILE *out,t_pull *pull,double t)
109 fprintf(out, "%.4f", t);
111 for(g=1; g<1+pull->ngrp; g++)
113 if (pull->eGeom == epullgPOS)
119 fprintf(out,"\t%g",pull->grp[g].f[d]);
125 fprintf(out,"\t%g",pull->grp[g].f_scal);
131 void pull_print_output(t_pull *pull, gmx_large_int_t step, double time)
133 if ((pull->nstxout != 0) && (step % pull->nstxout == 0))
135 pull_print_x(pull->out_x,pull,time);
138 if ((pull->nstfout != 0) && (step % pull->nstfout == 0))
140 pull_print_f(pull->out_f,pull,time);
144 static FILE *open_pull_out(const char *fn,t_pull *pull,const output_env_t oenv,
145 gmx_bool bCoord, unsigned long Flags)
149 char **setname,buf[10];
151 if(Flags & MD_APPENDFILES)
153 fp = gmx_fio_fopen(fn,"a+");
157 fp = gmx_fio_fopen(fn,"w+");
160 xvgr_header(fp,"Pull COM", "Time (ps)","Position (nm)",
165 xvgr_header(fp,"Pull force","Time (ps)","Force (kJ/mol/nm)",
169 snew(setname,(1+pull->ngrp)*DIM);
171 for(g=0; g<1+pull->ngrp; g++)
173 if (pull->grp[g].nat > 0 &&
174 (g > 0 || (bCoord && !PULL_CYL(pull))))
176 if (bCoord || pull->eGeom == epullgPOS)
184 sprintf(buf,"%d %s%c",g,"c",'X'+m);
185 setname[nsets] = strdup(buf);
194 sprintf(buf,"%d %s%c",
195 g,(bCoord && g > 0)?"d":"",'X'+m);
196 setname[nsets] = strdup(buf);
204 setname[nsets] = strdup(buf);
209 if (bCoord || nsets > 1)
211 xvgr_legend(fp,nsets,(const char**)setname,oenv);
213 for(g=0; g<nsets; g++)
223 /* Apply forces in a mass weighted fashion */
224 static void apply_forces_grp(t_pullgrp *pgrp, t_mdatoms * md,
226 dvec f_pull, int sign, rvec *f)
228 int i,ii,m,start,end;
232 end = md->homenr + start;
234 inv_wm = pgrp->wscale*pgrp->invtm;
236 for(i=0; i<pgrp->nat_loc; i++)
238 ii = pgrp->ind_loc[i];
239 wmass = md->massT[ii];
240 if (pgrp->weight_loc)
242 wmass *= pgrp->weight_loc[i];
247 f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
252 /* Apply forces in a mass weighted fashion */
253 static void apply_forces(t_pull * pull, t_mdatoms * md, gmx_ga2la_t ga2la,
259 for(i=1; i<pull->ngrp+1; i++)
261 pgrp = &(pull->grp[i]);
262 apply_forces_grp(pgrp,md,ga2la,pgrp->f,1,f);
263 if (pull->grp[0].nat)
267 apply_forces_grp(&(pull->dyna[i]),md,ga2la,pgrp->f,-1,f);
271 apply_forces_grp(&(pull->grp[0]),md,ga2la,pgrp->f,-1,f);
277 static double max_pull_distance2(const t_pull *pull,const t_pbc *pbc)
282 max_d2 = GMX_DOUBLE_MAX;
284 if (pull->eGeom != epullgDIRPBC)
286 for(m=0; m<pbc->ndim_ePBC; m++)
288 if (pull->dim[m] != 0)
290 max_d2 = min(max_d2,norm2(pbc->box[m]));
298 static void get_pullgrps_dr(const t_pull *pull,const t_pbc *pbc,int g,double t,
299 dvec xg,dvec xref,double max_dist2,
302 t_pullgrp *pref,*pgrp;
304 dvec xrefr,dref={0,0,0};
307 pgrp = &pull->grp[g];
309 copy_dvec(xref,xrefr);
311 if (pull->eGeom == epullgDIRPBC)
315 dref[m] = (pgrp->init[0] + pgrp->rate*t)*pull->grp[g].vec[m];
317 /* Add the reference position, so we use the correct periodic image */
318 dvec_inc(xrefr,dref);
321 pbc_dx_d(pbc, xg, xrefr, dr);
325 dr[m] *= pull->dim[m];
328 if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
330 gmx_fatal(FARGS,"Distance of pull group %d (%f nm) is larger than 0.49 times the box size (%f)",g,sqrt(dr2),sqrt(max_dist2));
333 if (pull->eGeom == epullgDIRPBC)
339 static void get_pullgrp_dr(const t_pull *pull,const t_pbc *pbc,int g,double t,
344 if (pull->eGeom == epullgDIRPBC)
350 md2 = max_pull_distance2(pull,pbc);
353 get_pullgrps_dr(pull,pbc,g,t,
355 PULL_CYL(pull) ? pull->dyna[g].x : pull->grp[0].x,
360 void get_pullgrp_distance(t_pull *pull,t_pbc *pbc,int g,double t,
363 static gmx_bool bWarned=FALSE; /* TODO: this should be fixed for thread-safety,
364 but is fairly benign */
370 pgrp = &pull->grp[g];
372 get_pullgrp_dr(pull,pbc,g,t,dr);
374 if (pull->eGeom == epullgPOS)
378 ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m];
383 ref[0] = pgrp->init[0] + pgrp->rate*t;
389 /* Pull along the vector between the com's */
390 if (ref[0] < 0 && !bWarned)
392 fprintf(stderr,"\nPull reference distance for group %d is negative (%f)\n",g,ref[0]);
398 /* With no vector we can not determine the direction for the force,
399 * so we set the force to zero.
405 /* Determine the deviation */
406 dev[0] = drs - ref[0];
416 inpr += pgrp->vec[m]*dr[m];
418 dev[0] = inpr - ref[0];
421 /* Determine the difference of dr and ref along each dimension */
424 dev[m] = (dr[m] - ref[m])*pull->dim[m];
430 void clear_pull_forces(t_pull *pull)
434 /* Zeroing the forces is only required for constraint pulling.
435 * It can happen that multiple constraint steps need to be applied
436 * and therefore the constraint forces need to be accumulated.
438 for(i=0; i<1+pull->ngrp; i++)
440 clear_dvec(pull->grp[i].f);
441 pull->grp[i].f_scal = 0;
445 /* Apply constraint using SHAKE */
446 static void do_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
448 gmx_bool bMaster, tensor vir,
452 dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
453 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
455 dvec *rinew; /* current 'new' position of group i */
456 dvec *rjnew; /* current 'new' position of group j */
459 double lambda, rm, mass, invdt=0;
460 gmx_bool bConverged_all,bConverged=FALSE;
461 int niter=0,g,ii,j,m,max_iter=100;
462 double q,a,b,c; /* for solving the quadratic equation,
463 see Num. Recipes in C ed 2 p. 184 */
464 dvec *dr; /* correction for group i */
465 dvec ref_dr; /* correction for group j */
466 dvec f; /* the pull force */
468 t_pullgrp *pdyna,*pgrp,*pref;
470 snew(r_ij,pull->ngrp+1);
473 snew(rjnew,pull->ngrp+1);
479 snew(dr,pull->ngrp+1);
480 snew(rinew,pull->ngrp+1);
482 /* copy the current unconstrained positions for use in iterations. We
483 iterate until rinew[i] and rjnew[j] obey the constraints. Then
484 rinew - pull.x_unc[i] is the correction dr to group i */
485 for(g=1; g<1+pull->ngrp; g++)
487 copy_dvec(pull->grp[g].xp,rinew[g]);
491 for(g=1; g<1+pull->ngrp; g++)
493 copy_dvec(pull->dyna[g].xp,rjnew[g]);
498 copy_dvec(pull->grp[0].xp,rjnew[0]);
501 /* Determine the constraint directions from the old positions */
502 for(g=1; g<1+pull->ngrp; g++)
504 get_pullgrp_dr(pull,pbc,g,t,r_ij[g]);
505 /* Store the difference vector at time t for printing */
506 copy_dvec(r_ij[g],pull->grp[g].dr);
509 fprintf(debug,"Pull group %d dr %f %f %f\n",
510 g,r_ij[g][XX],r_ij[g][YY],r_ij[g][ZZ]);
513 if (pull->eGeom == epullgDIR || pull->eGeom == epullgDIRPBC)
515 /* Select the component along vec */
519 a += pull->grp[g].vec[m]*r_ij[g][m];
523 r_ij[g][m] = a*pull->grp[g].vec[m];
528 bConverged_all = FALSE;
529 while (!bConverged_all && niter < max_iter)
531 bConverged_all = TRUE;
533 /* loop over all constraints */
534 for(g=1; g<1+pull->ngrp; g++)
536 pgrp = &pull->grp[g];
538 pref = &pull->dyna[g];
540 pref = &pull->grp[0];
542 /* Get the current difference vector */
543 get_pullgrps_dr(pull,pbc,g,t,rinew[g],rjnew[PULL_CYL(pull) ? g : 0],
546 if (pull->eGeom == epullgPOS)
550 ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m];
555 ref[0] = pgrp->init[0] + pgrp->rate*t;
556 /* Keep the compiler happy */
563 fprintf(debug,"Pull group %d, iteration %d\n",g,niter);
566 rm = 1.0/(pull->grp[g].invtm + pref->invtm);
573 gmx_fatal(FARGS,"The pull constraint reference distance for group %d is <= 0 (%f)",g,ref[0]);
576 a = diprod(r_ij[g],r_ij[g]);
577 b = diprod(unc_ij,r_ij[g])*2;
578 c = diprod(unc_ij,unc_ij) - dsqr(ref[0]);
582 q = -0.5*(b - sqrt(b*b - 4*a*c));
587 q = -0.5*(b + sqrt(b*b - 4*a*c));
594 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
598 /* The position corrections dr due to the constraints */
599 dsvmul(-lambda*rm*pgrp->invtm, r_ij[g], dr[g]);
600 dsvmul( lambda*rm*pref->invtm, r_ij[g], ref_dr);
605 /* A 1-dimensional constraint along a vector */
609 vec[m] = pgrp->vec[m];
610 a += unc_ij[m]*vec[m];
612 /* Select only the component along the vector */
613 dsvmul(a,vec,unc_ij);
617 fprintf(debug,"Pull inpr %e lambda: %e\n",a,lambda);
620 /* The position corrections dr due to the constraints */
621 dsvmul(-lambda*rm*pull->grp[g].invtm, vec, dr[g]);
622 dsvmul( lambda*rm* pref->invtm, vec,ref_dr);
629 lambda = r_ij[g][m] - ref[m];
630 /* The position corrections dr due to the constraints */
631 dr[g][m] = -lambda*rm*pull->grp[g].invtm;
632 ref_dr[m] = lambda*rm*pref->invtm;
646 j = (PULL_CYL(pull) ? g : 0);
647 get_pullgrps_dr(pull,pbc,g,t,rinew[g],rjnew[j],-1,tmp);
648 get_pullgrps_dr(pull,pbc,g,t,dr[g] ,ref_dr ,-1,tmp3);
650 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
651 rinew[g][0],rinew[g][1],rinew[g][2],
652 rjnew[j][0],rjnew[j][1],rjnew[j][2], dnorm(tmp));
653 if (pull->eGeom == epullgPOS)
656 "Pull ref %8.5f %8.5f %8.5f\n",
657 pgrp->vec[0],pgrp->vec[1],pgrp->vec[2]);
662 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f %8.5f %8.5f\n",
663 "","","","","","",ref[0],ref[1],ref[2]);
666 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
667 dr[g][0],dr[g][1],dr[g][2],
668 ref_dr[0],ref_dr[1],ref_dr[2],
671 "Pull cor %10.7f %10.7f %10.7f\n",
672 dr[g][0],dr[g][1],dr[g][2]);
675 /* Update the COMs with dr */
676 dvec_inc(rinew[g], dr[g]);
677 dvec_inc(rjnew[PULL_CYL(pull) ? g : 0],ref_dr);
680 /* Check if all constraints are fullfilled now */
681 for(g=1; g<1+pull->ngrp; g++)
683 pgrp = &pull->grp[g];
685 get_pullgrps_dr(pull,pbc,g,t,rinew[g],rjnew[PULL_CYL(pull) ? g : 0],
691 bConverged = fabs(dnorm(unc_ij) - ref[0]) < pull->constr_tol;
698 vec[m] = pgrp->vec[m];
700 inpr = diprod(unc_ij,vec);
701 dsvmul(inpr,vec,unc_ij);
703 fabs(diprod(unc_ij,vec) - ref[0]) < pull->constr_tol;
710 fabs(unc_ij[m] - ref[m]) >= pull->constr_tol)
722 fprintf(debug,"NOT CONVERGED YET: Group %d:"
723 "d_ref = %f %f %f, current d = %f\n",
724 g,ref[0],ref[1],ref[2],dnorm(unc_ij));
727 bConverged_all = FALSE;
732 /* if after all constraints are dealt with and bConverged is still TRUE
733 we're finished, if not we do another iteration */
735 if (niter > max_iter)
737 gmx_fatal(FARGS,"Too many iterations for constraint run: %d",niter);
740 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
747 /* update the normal groups */
748 for(g=1; g<1+pull->ngrp; g++)
750 pgrp = &pull->grp[g];
751 /* get the final dr and constraint force for group i */
752 dvec_sub(rinew[g],pgrp->xp,dr[g]);
753 /* select components of dr */
756 dr[g][m] *= pull->dim[m];
758 dsvmul(1.0/(pgrp->invtm*dt*dt),dr[g],f);
765 pgrp->f_scal += r_ij[g][m]*f[m]/dnorm(r_ij[g]);
773 pgrp->f_scal += pgrp->vec[m]*f[m];
780 if (vir && bMaster) {
781 /* Add the pull contribution to the virial */
786 vir[j][m] -= 0.5*f[j]*r_ij[g][m];
791 /* update the atom positions */
792 copy_dvec(dr[g],tmp);
793 for(j=0;j<pgrp->nat_loc;j++)
795 ii = pgrp->ind_loc[j];
796 if (pgrp->weight_loc)
798 dsvmul(pgrp->wscale*pgrp->weight_loc[j],dr[g],tmp);
808 v[ii][m] += invdt*tmp[m];
814 /* update the reference groups */
817 /* update the dynamic reference groups */
818 for(g=1; g<1+pull->ngrp; g++)
820 pdyna = &pull->dyna[g];
821 dvec_sub(rjnew[g],pdyna->xp,ref_dr);
822 /* select components of ref_dr */
825 ref_dr[m] *= pull->dim[m];
828 for(j=0;j<pdyna->nat_loc;j++)
830 /* reset the atoms with dr, weighted by w_i */
831 dsvmul(pdyna->wscale*pdyna->weight_loc[j],ref_dr,tmp);
832 ii = pdyna->ind_loc[j];
841 v[ii][m] += invdt*tmp[m];
849 pgrp = &pull->grp[0];
850 /* update the reference group */
851 dvec_sub(rjnew[0],pgrp->xp, ref_dr);
852 /* select components of ref_dr */
855 ref_dr[m] *= pull->dim[m];
858 copy_dvec(ref_dr,tmp);
859 for(j=0; j<pgrp->nat_loc;j++)
861 ii = pgrp->ind_loc[j];
862 if (pgrp->weight_loc)
864 dsvmul(pgrp->wscale*pgrp->weight_loc[j],ref_dr,tmp);
874 v[ii][m] += invdt*tmp[m];
880 /* finished! I hope. Give back some memory */
887 /* Pulling with a harmonic umbrella potential or constant force */
888 static void do_pull_pot(int ePull,
889 t_pull *pull, t_pbc *pbc, double t, real lambda,
890 real *V, tensor vir, real *dVdl)
898 /* loop over the groups that are being pulled */
901 for(g=1; g<1+pull->ngrp; g++)
903 pgrp = &pull->grp[g];
904 get_pullgrp_distance(pull,pbc,g,t,pgrp->dr,dev);
906 k = (1.0 - lambda)*pgrp->k + lambda*pgrp->kB;
907 dkdl = pgrp->kB - pgrp->k;
912 ndr = dnorm(pgrp->dr);
914 if (ePull == epullUMBRELLA)
916 pgrp->f_scal = -k*dev[0];
917 *V += 0.5* k*dsqr(dev[0]);
918 *dVdl += 0.5*dkdl*dsqr(dev[0]);
928 pgrp->f[m] = pgrp->f_scal*pgrp->dr[m]*invdr;
934 if (ePull == epullUMBRELLA)
936 pgrp->f_scal = -k*dev[0];
937 *V += 0.5* k*dsqr(dev[0]);
938 *dVdl += 0.5*dkdl*dsqr(dev[0]);
945 ndr += pgrp->vec[m]*pgrp->dr[m];
953 pgrp->f[m] = pgrp->f_scal*pgrp->vec[m];
959 if (ePull == epullUMBRELLA)
961 pgrp->f[m] = -k*dev[m];
962 *V += 0.5* k*dsqr(dev[m]);
963 *dVdl += 0.5*dkdl*dsqr(dev[m]);
967 pgrp->f[m] = -k*pull->dim[m];
968 *V += k*pgrp->dr[m]*pull->dim[m];
969 *dVdl += dkdl*pgrp->dr[m]*pull->dim[m];
977 /* Add the pull contribution to the virial */
982 vir[j][m] -= 0.5*pgrp->f[j]*pgrp->dr[m];
989 real pull_potential(int ePull,t_pull *pull, t_mdatoms *md, t_pbc *pbc,
990 t_commrec *cr, double t, real lambda,
991 rvec *x, rvec *f, tensor vir, real *dvdlambda)
995 pull_calc_coms(cr,pull,md,pbc,t,x,NULL);
997 do_pull_pot(ePull,pull,pbc,t,lambda,
998 &V,pull->bVirial && MASTER(cr) ? vir : NULL,&dVdl);
1000 /* Distribute forces over pulled groups */
1001 apply_forces(pull, md, DOMAINDECOMP(cr) ? cr->dd->ga2la : NULL, f);
1007 return (MASTER(cr) ? V : 0.0);
1010 void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1011 t_commrec *cr, double dt, double t,
1012 rvec *x, rvec *xp, rvec *v, tensor vir)
1014 pull_calc_coms(cr,pull,md,pbc,t,x,xp);
1016 do_constraint(pull,md,pbc,xp,v,pull->bVirial && MASTER(cr),vir,dt,t);
1019 static void make_local_pull_group(gmx_ga2la_t ga2la,
1020 t_pullgrp *pg,int start,int end)
1025 for(i=0; i<pg->nat; i++) {
1028 if (!ga2la_get_home(ga2la,ii,&ii)) {
1032 if (ii >= start && ii < end) {
1033 /* This is a home atom, add it to the local pull group */
1034 if (pg->nat_loc >= pg->nalloc_loc) {
1035 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1036 srenew(pg->ind_loc,pg->nalloc_loc);
1037 if (pg->epgrppbc == epgrppbcCOS || pg->weight) {
1038 srenew(pg->weight_loc,pg->nalloc_loc);
1041 pg->ind_loc[pg->nat_loc] = ii;
1043 pg->weight_loc[pg->nat_loc] = pg->weight[i];
1050 void dd_make_local_pull_groups(gmx_domdec_t *dd,t_pull *pull,t_mdatoms *md)
1061 if (pull->grp[0].nat > 0)
1062 make_local_pull_group(ga2la,&pull->grp[0],md->start,md->start+md->homenr);
1063 for(g=1; g<1+pull->ngrp; g++)
1064 make_local_pull_group(ga2la,&pull->grp[g],md->start,md->start+md->homenr);
1067 static void init_pull_group_index(FILE *fplog,t_commrec *cr,
1069 int g,t_pullgrp *pg,ivec pulldims,
1070 gmx_mtop_t *mtop,t_inputrec *ir, real lambda)
1072 int i,ii,d,nfrozen,ndim;
1074 double tmass,wmass,wwmass;
1076 gmx_ga2la_t ga2la=NULL;
1077 gmx_groups_t *groups;
1078 gmx_mtop_atomlookup_t alook;
1081 bDomDec = (cr && DOMAINDECOMP(cr));
1083 ga2la = cr->dd->ga2la;
1086 if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD) {
1087 /* There are no masses in the integrator.
1088 * But we still want to have the correct mass-weighted COMs.
1089 * So we store the real masses in the weights.
1090 * We do not set nweight, so these weights do not end up in the tpx file.
1092 if (pg->nweight == 0) {
1093 snew(pg->weight,pg->nat);
1097 if (cr && PAR(cr)) {
1101 pg->weight_loc = NULL;
1103 pg->nat_loc = pg->nat;
1104 pg->ind_loc = pg->ind;
1105 if (pg->epgrppbc == epgrppbcCOS) {
1106 snew(pg->weight_loc,pg->nat);
1108 pg->weight_loc = pg->weight;
1112 groups = &mtop->groups;
1114 alook = gmx_mtop_atomlookup_init(mtop);
1120 for(i=0; i<pg->nat; i++) {
1122 gmx_mtop_atomnr_to_atom(alook,ii,&atom);
1123 if (cr && PAR(cr) && !bDomDec && ii >= start && ii < end)
1124 pg->ind_loc[pg->nat_loc++] = ii;
1125 if (ir->opts.nFreeze) {
1126 for(d=0; d<DIM; d++)
1127 if (pulldims[d] && ir->opts.nFreeze[ggrpnr(groups,egcFREEZE,ii)][d])
1130 if (ir->efep == efepNO) {
1133 m = (1 - lambda)*atom->m + lambda*atom->mB;
1135 if (pg->nweight > 0) {
1140 if (EI_ENERGY_MINIMIZATION(ir->eI)) {
1141 /* Move the mass to the weight */
1145 } else if (ir->eI == eiBD) {
1147 mbd = ir->bd_fric*ir->delta_t;
1149 if (groups->grpnr[egcTC] == NULL) {
1150 mbd = ir->delta_t/ir->opts.tau_t[0];
1152 mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1164 gmx_mtop_atomlookup_destroy(alook);
1167 gmx_fatal(FARGS,"The total%s mass of pull group %d is zero",
1168 pg->weight ? " weighted" : "",g);
1172 "Pull group %d: %5d atoms, mass %9.3f",g,pg->nat,tmass);
1173 if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD) {
1174 fprintf(fplog,", weighted mass %9.3f",wmass*wmass/wwmass);
1176 if (pg->epgrppbc == epgrppbcCOS) {
1177 fprintf(fplog,", cosine weighting will be used");
1179 fprintf(fplog,"\n");
1183 /* A value > 0 signals not frozen, it is updated later */
1187 for(d=0; d<DIM; d++)
1188 ndim += pulldims[d]*pg->nat;
1189 if (fplog && nfrozen > 0 && nfrozen < ndim) {
1191 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1192 " that are subject to pulling are frozen.\n"
1193 " For pulling the whole group will be frozen.\n\n",
1201 void init_pull(FILE *fplog,t_inputrec *ir,int nfile,const t_filenm fnm[],
1202 gmx_mtop_t *mtop,t_commrec *cr,const output_env_t oenv, real lambda,
1203 gmx_bool bOutFile, unsigned long Flags)
1207 int g,start=0,end=0,m;
1212 pull->ePBC = ir->ePBC;
1215 case epbcNONE: pull->npbcdim = 0; break;
1216 case epbcXY: pull->npbcdim = 2; break;
1217 default: pull->npbcdim = 3; break;
1222 fprintf(fplog,"\nWill apply %s COM pulling in geometry '%s'\n",
1223 EPULLTYPE(ir->ePull),EPULLGEOM(pull->eGeom));
1224 if (pull->grp[0].nat > 0)
1226 fprintf(fplog,"between a reference group and %d group%s\n",
1227 pull->ngrp,pull->ngrp==1 ? "" : "s");
1231 fprintf(fplog,"with an absolute reference on %d group%s\n",
1232 pull->ngrp,pull->ngrp==1 ? "" : "s");
1235 for(g=0; g<pull->ngrp+1; g++)
1237 if (pull->grp[g].nat > 1 &&
1238 pull->grp[g].pbcatom < 0)
1240 /* We are using cosine weighting */
1241 fprintf(fplog,"Cosine weighting is used for group %d\n",g);
1247 please_cite(fplog,"Engin2010");
1251 /* We always add the virial contribution,
1252 * except for geometry = direction_periodic where this is impossible.
1254 pull->bVirial = (pull->eGeom != epullgDIRPBC);
1255 if (getenv("GMX_NO_PULLVIR") != NULL)
1259 fprintf(fplog,"Found env. var., will not add the virial contribution of the COM pull forces\n");
1261 pull->bVirial = FALSE;
1264 if (cr && PARTDECOMP(cr))
1266 pd_at_range(cr,&start,&end);
1270 pull->dbuf_cyl=NULL;
1271 pull->bRefAt = FALSE;
1273 for(g=0; g<pull->ngrp+1; g++)
1275 pgrp = &pull->grp[g];
1276 pgrp->epgrppbc = epgrppbcNONE;
1279 /* Determine if we need to take PBC into account for calculating
1280 * the COM's of the pull groups.
1282 for(m=0; m<pull->npbcdim; m++)
1284 if (pull->dim[m] && pgrp->nat > 1)
1286 if (pgrp->pbcatom >= 0)
1288 pgrp->epgrppbc = epgrppbcREFAT;
1289 pull->bRefAt = TRUE;
1295 gmx_fatal(FARGS,"Pull groups can not have relative weights and cosine weighting at same time");
1297 pgrp->epgrppbc = epgrppbcCOS;
1298 if (pull->cosdim >= 0 && pull->cosdim != m)
1300 gmx_fatal(FARGS,"Can only use cosine weighting with pulling in one dimension (use mdp option pull_dim)");
1306 /* Set the indices */
1307 init_pull_group_index(fplog,cr,start,end,g,pgrp,pull->dim,mtop,ir,lambda);
1308 if (PULL_CYL(pull) && pgrp->invtm == 0)
1310 gmx_fatal(FARGS,"Can not have frozen atoms in a cylinder pull group");
1315 /* Absolute reference, set the inverse mass to zero */
1321 /* if we use dynamic reference groups, do some initialising for them */
1324 if (pull->grp[0].nat == 0)
1326 gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n");
1328 snew(pull->dyna,pull->ngrp+1);
1331 /* Only do I/O when we are doing dynamics and if we are the MASTER */
1336 if (pull->nstxout > 0)
1338 pull->out_x = open_pull_out(opt2fn("-px",nfile,fnm),pull,oenv,TRUE,Flags);
1340 if (pull->nstfout > 0)
1342 pull->out_f = open_pull_out(opt2fn("-pf",nfile,fnm),pull,oenv,
1348 void finish_pull(FILE *fplog,t_pull *pull)
1352 gmx_fio_fclose(pull->out_x);
1356 gmx_fio_fclose(pull->out_f);