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"
64 static void pull_print_x_grp(FILE *out
,gmx_bool bRef
,ivec dim
,t_pullgrp
*pgrp
)
72 fprintf(out
,"\t%g",bRef
? pgrp
->x
[m
] : pgrp
->dr
[m
]);
77 static void pull_print_x(FILE *out
,t_pull
*pull
,double t
)
81 fprintf(out
, "%.4f", t
);
85 for (g
=1; g
<1+pull
->ngrp
; g
++)
87 pull_print_x_grp(out
,TRUE
,pull
->dim
,&pull
->dyna
[g
]);
88 pull_print_x_grp(out
,FALSE
,pull
->dim
,&pull
->grp
[g
]);
93 for (g
=0; g
<1+pull
->ngrp
; g
++)
95 if (pull
->grp
[g
].nat
> 0)
97 pull_print_x_grp(out
,g
==0,pull
->dim
,&pull
->grp
[g
]);
104 static void pull_print_f(FILE *out
,t_pull
*pull
,double t
)
108 fprintf(out
, "%.4f", t
);
110 for(g
=1; g
<1+pull
->ngrp
; g
++)
112 if (pull
->eGeom
== epullgPOS
)
118 fprintf(out
,"\t%g",pull
->grp
[g
].f
[d
]);
124 fprintf(out
,"\t%g",pull
->grp
[g
].f_scal
);
130 void pull_print_output(t_pull
*pull
, gmx_large_int_t step
, double time
)
132 if ((pull
->nstxout
!= 0) && (step
% pull
->nstxout
== 0))
134 pull_print_x(pull
->out_x
,pull
,time
);
137 if ((pull
->nstfout
!= 0) && (step
% pull
->nstfout
== 0))
139 pull_print_f(pull
->out_f
,pull
,time
);
143 static FILE *open_pull_out(const char *fn
,t_pull
*pull
,const output_env_t oenv
,
144 gmx_bool bCoord
, unsigned long Flags
)
148 char **setname
,buf
[10];
150 if(Flags
& MD_APPENDFILES
)
152 fp
= gmx_fio_fopen(fn
,"a+");
156 fp
= gmx_fio_fopen(fn
,"w+");
159 xvgr_header(fp
,"Pull COM", "Time (ps)","Position (nm)",
164 xvgr_header(fp
,"Pull force","Time (ps)","Force (kJ/mol/nm)",
168 snew(setname
,(1+pull
->ngrp
)*DIM
);
170 for(g
=0; g
<1+pull
->ngrp
; g
++)
172 if (pull
->grp
[g
].nat
> 0 &&
173 (g
> 0 || (bCoord
&& !PULL_CYL(pull
))))
175 if (bCoord
|| pull
->eGeom
== epullgPOS
)
183 sprintf(buf
,"%d %s%c",g
,"c",'X'+m
);
184 setname
[nsets
] = strdup(buf
);
193 sprintf(buf
,"%d %s%c",
194 g
,(bCoord
&& g
> 0)?"d":"",'X'+m
);
195 setname
[nsets
] = strdup(buf
);
203 setname
[nsets
] = strdup(buf
);
208 if (bCoord
|| nsets
> 1)
210 xvgr_legend(fp
,nsets
,(const char**)setname
,oenv
);
212 for(g
=0; g
<nsets
; g
++)
222 /* Apply forces in a mass weighted fashion */
223 static void apply_forces_grp(t_pullgrp
*pgrp
, t_mdatoms
* md
,
225 dvec f_pull
, int sign
, rvec
*f
)
227 int i
,ii
,m
,start
,end
;
231 end
= md
->homenr
+ start
;
233 inv_wm
= pgrp
->wscale
*pgrp
->invtm
;
235 for(i
=0; i
<pgrp
->nat_loc
; i
++)
237 ii
= pgrp
->ind_loc
[i
];
238 wmass
= md
->massT
[ii
];
239 if (pgrp
->weight_loc
)
241 wmass
*= pgrp
->weight_loc
[i
];
246 f
[ii
][m
] += sign
* wmass
* f_pull
[m
] * inv_wm
;
251 /* Apply forces in a mass weighted fashion */
252 static void apply_forces(t_pull
* pull
, t_mdatoms
* md
, gmx_ga2la_t ga2la
,
258 for(i
=1; i
<pull
->ngrp
+1; i
++)
260 pgrp
= &(pull
->grp
[i
]);
261 apply_forces_grp(pgrp
,md
,ga2la
,pgrp
->f
,1,f
);
262 if (pull
->grp
[0].nat
)
266 apply_forces_grp(&(pull
->dyna
[i
]),md
,ga2la
,pgrp
->f
,-1,f
);
270 apply_forces_grp(&(pull
->grp
[0]),md
,ga2la
,pgrp
->f
,-1,f
);
276 static double max_pull_distance2(const t_pull
*pull
,const t_pbc
*pbc
)
281 max_d2
= GMX_DOUBLE_MAX
;
283 if (pull
->eGeom
!= epullgDIRPBC
)
285 for(m
=0; m
<pbc
->ndim_ePBC
; m
++)
287 if (pull
->dim
[m
] != 0)
289 max_d2
= min(max_d2
,norm2(pbc
->box
[m
]));
297 static void get_pullgrps_dr(const t_pull
*pull
,const t_pbc
*pbc
,int g
,double t
,
298 dvec xg
,dvec xref
,double max_dist2
,
301 t_pullgrp
*pref
,*pgrp
;
303 dvec xrefr
,dref
={0,0,0};
306 pgrp
= &pull
->grp
[g
];
308 copy_dvec(xref
,xrefr
);
310 if (pull
->eGeom
== epullgDIRPBC
)
314 dref
[m
] = (pgrp
->init
[0] + pgrp
->rate
*t
)*pull
->grp
[g
].vec
[m
];
316 /* Add the reference position, so we use the correct periodic image */
317 dvec_inc(xrefr
,dref
);
320 pbc_dx_d(pbc
, xg
, xrefr
, dr
);
324 dr
[m
] *= pull
->dim
[m
];
327 if (max_dist2
>= 0 && dr2
> 0.98*0.98*max_dist2
)
329 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
));
332 if (pull
->eGeom
== epullgDIRPBC
)
338 static void get_pullgrp_dr(const t_pull
*pull
,const t_pbc
*pbc
,int g
,double t
,
343 if (pull
->eGeom
== epullgDIRPBC
)
349 md2
= max_pull_distance2(pull
,pbc
);
352 get_pullgrps_dr(pull
,pbc
,g
,t
,
354 PULL_CYL(pull
) ? pull
->dyna
[g
].x
: pull
->grp
[0].x
,
359 void get_pullgrp_distance(t_pull
*pull
,t_pbc
*pbc
,int g
,double t
,
362 static gmx_bool bWarned
=FALSE
; /* TODO: this should be fixed for thread-safety,
363 but is fairly benign */
369 pgrp
= &pull
->grp
[g
];
371 get_pullgrp_dr(pull
,pbc
,g
,t
,dr
);
373 if (pull
->eGeom
== epullgPOS
)
377 ref
[m
] = pgrp
->init
[m
] + pgrp
->rate
*t
*pgrp
->vec
[m
];
382 ref
[0] = pgrp
->init
[0] + pgrp
->rate
*t
;
388 /* Pull along the vector between the com's */
389 if (ref
[0] < 0 && !bWarned
)
391 fprintf(stderr
,"\nPull reference distance for group %d is negative (%f)\n",g
,ref
[0]);
397 /* With no vector we can not determine the direction for the force,
398 * so we set the force to zero.
404 /* Determine the deviation */
405 dev
[0] = drs
- ref
[0];
415 inpr
+= pgrp
->vec
[m
]*dr
[m
];
417 dev
[0] = inpr
- ref
[0];
420 /* Determine the difference of dr and ref along each dimension */
423 dev
[m
] = (dr
[m
] - ref
[m
])*pull
->dim
[m
];
429 void clear_pull_forces(t_pull
*pull
)
433 /* Zeroing the forces is only required for constraint pulling.
434 * It can happen that multiple constraint steps need to be applied
435 * and therefore the constraint forces need to be accumulated.
437 for(i
=0; i
<1+pull
->ngrp
; i
++)
439 clear_dvec(pull
->grp
[i
].f
);
440 pull
->grp
[i
].f_scal
= 0;
444 /* Apply constraint using SHAKE */
445 static void do_constraint(t_pull
*pull
, t_mdatoms
*md
, t_pbc
*pbc
,
447 gmx_bool bMaster
, tensor vir
,
451 dvec
*r_ij
; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
452 dvec unc_ij
; /* xp[i] com of i this step, before constr. -> unc_ij */
454 dvec
*rinew
; /* current 'new' position of group i */
455 dvec
*rjnew
; /* current 'new' position of group j */
458 double lambda
, rm
, mass
, invdt
=0;
459 gmx_bool bConverged_all
,bConverged
=FALSE
;
460 int niter
=0,g
,ii
,j
,m
,max_iter
=100;
461 double q
,a
,b
,c
; /* for solving the quadratic equation,
462 see Num. Recipes in C ed 2 p. 184 */
463 dvec
*dr
; /* correction for group i */
464 dvec ref_dr
; /* correction for group j */
465 dvec f
; /* the pull force */
467 t_pullgrp
*pdyna
,*pgrp
,*pref
;
469 snew(r_ij
,pull
->ngrp
+1);
472 snew(rjnew
,pull
->ngrp
+1);
478 snew(dr
,pull
->ngrp
+1);
479 snew(rinew
,pull
->ngrp
+1);
481 /* copy the current unconstrained positions for use in iterations. We
482 iterate until rinew[i] and rjnew[j] obey the constraints. Then
483 rinew - pull.x_unc[i] is the correction dr to group i */
484 for(g
=1; g
<1+pull
->ngrp
; g
++)
486 copy_dvec(pull
->grp
[g
].xp
,rinew
[g
]);
490 for(g
=1; g
<1+pull
->ngrp
; g
++)
492 copy_dvec(pull
->dyna
[g
].xp
,rjnew
[g
]);
497 copy_dvec(pull
->grp
[0].xp
,rjnew
[0]);
500 /* Determine the constraint directions from the old positions */
501 for(g
=1; g
<1+pull
->ngrp
; g
++)
503 get_pullgrp_dr(pull
,pbc
,g
,t
,r_ij
[g
]);
504 /* Store the difference vector at time t for printing */
505 copy_dvec(r_ij
[g
],pull
->grp
[g
].dr
);
508 fprintf(debug
,"Pull group %d dr %f %f %f\n",
509 g
,r_ij
[g
][XX
],r_ij
[g
][YY
],r_ij
[g
][ZZ
]);
512 if (pull
->eGeom
== epullgDIR
|| pull
->eGeom
== epullgDIRPBC
)
514 /* Select the component along vec */
518 a
+= pull
->grp
[g
].vec
[m
]*r_ij
[g
][m
];
522 r_ij
[g
][m
] = a
*pull
->grp
[g
].vec
[m
];
527 bConverged_all
= FALSE
;
528 while (!bConverged_all
&& niter
< max_iter
)
530 bConverged_all
= TRUE
;
532 /* loop over all constraints */
533 for(g
=1; g
<1+pull
->ngrp
; g
++)
535 pgrp
= &pull
->grp
[g
];
537 pref
= &pull
->dyna
[g
];
539 pref
= &pull
->grp
[0];
541 /* Get the current difference vector */
542 get_pullgrps_dr(pull
,pbc
,g
,t
,rinew
[g
],rjnew
[PULL_CYL(pull
) ? g
: 0],
545 if (pull
->eGeom
== epullgPOS
)
549 ref
[m
] = pgrp
->init
[m
] + pgrp
->rate
*t
*pgrp
->vec
[m
];
554 ref
[0] = pgrp
->init
[0] + pgrp
->rate
*t
;
555 /* Keep the compiler happy */
562 fprintf(debug
,"Pull group %d, iteration %d\n",g
,niter
);
565 rm
= 1.0/(pull
->grp
[g
].invtm
+ pref
->invtm
);
572 gmx_fatal(FARGS
,"The pull constraint reference distance for group %d is <= 0 (%f)",g
,ref
[0]);
575 a
= diprod(r_ij
[g
],r_ij
[g
]);
576 b
= diprod(unc_ij
,r_ij
[g
])*2;
577 c
= diprod(unc_ij
,unc_ij
) - dsqr(ref
[0]);
581 q
= -0.5*(b
- sqrt(b
*b
- 4*a
*c
));
586 q
= -0.5*(b
+ sqrt(b
*b
- 4*a
*c
));
593 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
597 /* The position corrections dr due to the constraints */
598 dsvmul(-lambda
*rm
*pgrp
->invtm
, r_ij
[g
], dr
[g
]);
599 dsvmul( lambda
*rm
*pref
->invtm
, r_ij
[g
], ref_dr
);
604 /* A 1-dimensional constraint along a vector */
608 vec
[m
] = pgrp
->vec
[m
];
609 a
+= unc_ij
[m
]*vec
[m
];
611 /* Select only the component along the vector */
612 dsvmul(a
,vec
,unc_ij
);
616 fprintf(debug
,"Pull inpr %e lambda: %e\n",a
,lambda
);
619 /* The position corrections dr due to the constraints */
620 dsvmul(-lambda
*rm
*pull
->grp
[g
].invtm
, vec
, dr
[g
]);
621 dsvmul( lambda
*rm
* pref
->invtm
, vec
,ref_dr
);
628 lambda
= r_ij
[g
][m
] - ref
[m
];
629 /* The position corrections dr due to the constraints */
630 dr
[g
][m
] = -lambda
*rm
*pull
->grp
[g
].invtm
;
631 ref_dr
[m
] = lambda
*rm
*pref
->invtm
;
645 j
= (PULL_CYL(pull
) ? g
: 0);
646 get_pullgrps_dr(pull
,pbc
,g
,t
,rinew
[g
],rjnew
[j
],-1,tmp
);
647 get_pullgrps_dr(pull
,pbc
,g
,t
,dr
[g
] ,ref_dr
,-1,tmp3
);
649 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
650 rinew
[g
][0],rinew
[g
][1],rinew
[g
][2],
651 rjnew
[j
][0],rjnew
[j
][1],rjnew
[j
][2], dnorm(tmp
));
652 if (pull
->eGeom
== epullgPOS
)
655 "Pull ref %8.5f %8.5f %8.5f\n",
656 pgrp
->vec
[0],pgrp
->vec
[1],pgrp
->vec
[2]);
661 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f %8.5f %8.5f\n",
662 "","","","","","",ref
[0],ref
[1],ref
[2]);
665 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
666 dr
[g
][0],dr
[g
][1],dr
[g
][2],
667 ref_dr
[0],ref_dr
[1],ref_dr
[2],
670 "Pull cor %10.7f %10.7f %10.7f\n",
671 dr
[g
][0],dr
[g
][1],dr
[g
][2]);
674 /* Update the COMs with dr */
675 dvec_inc(rinew
[g
], dr
[g
]);
676 dvec_inc(rjnew
[PULL_CYL(pull
) ? g
: 0],ref_dr
);
679 /* Check if all constraints are fullfilled now */
680 for(g
=1; g
<1+pull
->ngrp
; g
++)
682 pgrp
= &pull
->grp
[g
];
684 get_pullgrps_dr(pull
,pbc
,g
,t
,rinew
[g
],rjnew
[PULL_CYL(pull
) ? g
: 0],
690 bConverged
= fabs(dnorm(unc_ij
) - ref
[0]) < pull
->constr_tol
;
697 vec
[m
] = pgrp
->vec
[m
];
699 inpr
= diprod(unc_ij
,vec
);
700 dsvmul(inpr
,vec
,unc_ij
);
702 fabs(diprod(unc_ij
,vec
) - ref
[0]) < pull
->constr_tol
;
709 fabs(unc_ij
[m
] - ref
[m
]) >= pull
->constr_tol
)
721 fprintf(debug
,"NOT CONVERGED YET: Group %d:"
722 "d_ref = %f %f %f, current d = %f\n",
723 g
,ref
[0],ref
[1],ref
[2],dnorm(unc_ij
));
726 bConverged_all
= FALSE
;
731 /* if after all constraints are dealt with and bConverged is still TRUE
732 we're finished, if not we do another iteration */
734 if (niter
> max_iter
)
736 gmx_fatal(FARGS
,"Too many iterations for constraint run: %d",niter
);
739 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
746 /* update the normal groups */
747 for(g
=1; g
<1+pull
->ngrp
; g
++)
749 pgrp
= &pull
->grp
[g
];
750 /* get the final dr and constraint force for group i */
751 dvec_sub(rinew
[g
],pgrp
->xp
,dr
[g
]);
752 /* select components of dr */
755 dr
[g
][m
] *= pull
->dim
[m
];
757 dsvmul(1.0/(pgrp
->invtm
*dt
*dt
),dr
[g
],f
);
764 pgrp
->f_scal
+= r_ij
[g
][m
]*f
[m
]/dnorm(r_ij
[g
]);
772 pgrp
->f_scal
+= pgrp
->vec
[m
]*f
[m
];
779 if (vir
&& bMaster
) {
780 /* Add the pull contribution to the virial */
785 vir
[j
][m
] -= 0.5*f
[j
]*r_ij
[g
][m
];
790 /* update the atom positions */
791 copy_dvec(dr
[g
],tmp
);
792 for(j
=0;j
<pgrp
->nat_loc
;j
++)
794 ii
= pgrp
->ind_loc
[j
];
795 if (pgrp
->weight_loc
)
797 dsvmul(pgrp
->wscale
*pgrp
->weight_loc
[j
],dr
[g
],tmp
);
807 v
[ii
][m
] += invdt
*tmp
[m
];
813 /* update the reference groups */
816 /* update the dynamic reference groups */
817 for(g
=1; g
<1+pull
->ngrp
; g
++)
819 pdyna
= &pull
->dyna
[g
];
820 dvec_sub(rjnew
[g
],pdyna
->xp
,ref_dr
);
821 /* select components of ref_dr */
824 ref_dr
[m
] *= pull
->dim
[m
];
827 for(j
=0;j
<pdyna
->nat_loc
;j
++)
829 /* reset the atoms with dr, weighted by w_i */
830 dsvmul(pdyna
->wscale
*pdyna
->weight_loc
[j
],ref_dr
,tmp
);
831 ii
= pdyna
->ind_loc
[j
];
840 v
[ii
][m
] += invdt
*tmp
[m
];
848 pgrp
= &pull
->grp
[0];
849 /* update the reference group */
850 dvec_sub(rjnew
[0],pgrp
->xp
, ref_dr
);
851 /* select components of ref_dr */
854 ref_dr
[m
] *= pull
->dim
[m
];
857 copy_dvec(ref_dr
,tmp
);
858 for(j
=0; j
<pgrp
->nat_loc
;j
++)
860 ii
= pgrp
->ind_loc
[j
];
861 if (pgrp
->weight_loc
)
863 dsvmul(pgrp
->wscale
*pgrp
->weight_loc
[j
],ref_dr
,tmp
);
873 v
[ii
][m
] += invdt
*tmp
[m
];
879 /* finished! I hope. Give back some memory */
886 /* Pulling with a harmonic umbrella potential or constant force */
887 static void do_pull_pot(int ePull
,
888 t_pull
*pull
, t_pbc
*pbc
, double t
, real lambda
,
889 real
*V
, tensor vir
, real
*dVdl
)
897 /* loop over the groups that are being pulled */
900 for(g
=1; g
<1+pull
->ngrp
; g
++)
902 pgrp
= &pull
->grp
[g
];
903 get_pullgrp_distance(pull
,pbc
,g
,t
,pgrp
->dr
,dev
);
905 k
= (1.0 - lambda
)*pgrp
->k
+ lambda
*pgrp
->kB
;
906 dkdl
= pgrp
->kB
- pgrp
->k
;
911 ndr
= dnorm(pgrp
->dr
);
913 if (ePull
== epullUMBRELLA
)
915 pgrp
->f_scal
= -k
*dev
[0];
916 *V
+= 0.5* k
*dsqr(dev
[0]);
917 *dVdl
+= 0.5*dkdl
*dsqr(dev
[0]);
927 pgrp
->f
[m
] = pgrp
->f_scal
*pgrp
->dr
[m
]*invdr
;
933 if (ePull
== epullUMBRELLA
)
935 pgrp
->f_scal
= -k
*dev
[0];
936 *V
+= 0.5* k
*dsqr(dev
[0]);
937 *dVdl
+= 0.5*dkdl
*dsqr(dev
[0]);
944 ndr
+= pgrp
->vec
[m
]*pgrp
->dr
[m
];
952 pgrp
->f
[m
] = pgrp
->f_scal
*pgrp
->vec
[m
];
958 if (ePull
== epullUMBRELLA
)
960 pgrp
->f
[m
] = -k
*dev
[m
];
961 *V
+= 0.5* k
*dsqr(dev
[m
]);
962 *dVdl
+= 0.5*dkdl
*dsqr(dev
[m
]);
966 pgrp
->f
[m
] = -k
*pull
->dim
[m
];
967 *V
+= k
*pgrp
->dr
[m
]*pull
->dim
[m
];
968 *dVdl
+= dkdl
*pgrp
->dr
[m
]*pull
->dim
[m
];
976 /* Add the pull contribution to the virial */
981 vir
[j
][m
] -= 0.5*pgrp
->f
[j
]*pgrp
->dr
[m
];
988 real
pull_potential(int ePull
,t_pull
*pull
, t_mdatoms
*md
, t_pbc
*pbc
,
989 t_commrec
*cr
, double t
, real lambda
,
990 rvec
*x
, rvec
*f
, tensor vir
, real
*dvdlambda
)
994 pull_calc_coms(cr
,pull
,md
,pbc
,t
,x
,NULL
);
996 do_pull_pot(ePull
,pull
,pbc
,t
,lambda
,
997 &V
,pull
->bVirial
&& MASTER(cr
) ? vir
: NULL
,&dVdl
);
999 /* Distribute forces over pulled groups */
1000 apply_forces(pull
, md
, DOMAINDECOMP(cr
) ? cr
->dd
->ga2la
: NULL
, f
);
1006 return (MASTER(cr
) ? V
: 0.0);
1009 void pull_constraint(t_pull
*pull
, t_mdatoms
*md
, t_pbc
*pbc
,
1010 t_commrec
*cr
, double dt
, double t
,
1011 rvec
*x
, rvec
*xp
, rvec
*v
, tensor vir
)
1013 pull_calc_coms(cr
,pull
,md
,pbc
,t
,x
,xp
);
1015 do_constraint(pull
,md
,pbc
,xp
,v
,pull
->bVirial
&& MASTER(cr
),vir
,dt
,t
);
1018 static void make_local_pull_group(gmx_ga2la_t ga2la
,
1019 t_pullgrp
*pg
,int start
,int end
)
1024 for(i
=0; i
<pg
->nat
; i
++) {
1027 if (!ga2la_get_home(ga2la
,ii
,&ii
)) {
1031 if (ii
>= start
&& ii
< end
) {
1032 /* This is a home atom, add it to the local pull group */
1033 if (pg
->nat_loc
>= pg
->nalloc_loc
) {
1034 pg
->nalloc_loc
= over_alloc_dd(pg
->nat_loc
+1);
1035 srenew(pg
->ind_loc
,pg
->nalloc_loc
);
1036 if (pg
->epgrppbc
== epgrppbcCOS
|| pg
->weight
) {
1037 srenew(pg
->weight_loc
,pg
->nalloc_loc
);
1040 pg
->ind_loc
[pg
->nat_loc
] = ii
;
1042 pg
->weight_loc
[pg
->nat_loc
] = pg
->weight
[i
];
1049 void dd_make_local_pull_groups(gmx_domdec_t
*dd
,t_pull
*pull
,t_mdatoms
*md
)
1060 if (pull
->grp
[0].nat
> 0)
1061 make_local_pull_group(ga2la
,&pull
->grp
[0],md
->start
,md
->start
+md
->homenr
);
1062 for(g
=1; g
<1+pull
->ngrp
; g
++)
1063 make_local_pull_group(ga2la
,&pull
->grp
[g
],md
->start
,md
->start
+md
->homenr
);
1066 static void init_pull_group_index(FILE *fplog
,t_commrec
*cr
,
1068 int g
,t_pullgrp
*pg
,ivec pulldims
,
1069 gmx_mtop_t
*mtop
,t_inputrec
*ir
, real lambda
)
1071 int i
,ii
,d
,nfrozen
,ndim
;
1073 double tmass
,wmass
,wwmass
;
1075 gmx_ga2la_t ga2la
=NULL
;
1076 gmx_groups_t
*groups
;
1077 gmx_mtop_atomlookup_t alook
;
1080 bDomDec
= (cr
&& DOMAINDECOMP(cr
));
1082 ga2la
= cr
->dd
->ga2la
;
1085 if (EI_ENERGY_MINIMIZATION(ir
->eI
) || ir
->eI
== eiBD
) {
1086 /* There are no masses in the integrator.
1087 * But we still want to have the correct mass-weighted COMs.
1088 * So we store the real masses in the weights.
1089 * We do not set nweight, so these weights do not end up in the tpx file.
1091 if (pg
->nweight
== 0) {
1092 snew(pg
->weight
,pg
->nat
);
1096 if (cr
&& PAR(cr
)) {
1100 pg
->weight_loc
= NULL
;
1102 pg
->nat_loc
= pg
->nat
;
1103 pg
->ind_loc
= pg
->ind
;
1104 if (pg
->epgrppbc
== epgrppbcCOS
) {
1105 snew(pg
->weight_loc
,pg
->nat
);
1107 pg
->weight_loc
= pg
->weight
;
1111 groups
= &mtop
->groups
;
1113 alook
= gmx_mtop_atomlookup_init(mtop
);
1119 for(i
=0; i
<pg
->nat
; i
++) {
1121 gmx_mtop_atomnr_to_atom(alook
,ii
,&atom
);
1122 if (cr
&& PAR(cr
) && !bDomDec
&& ii
>= start
&& ii
< end
)
1123 pg
->ind_loc
[pg
->nat_loc
++] = ii
;
1124 if (ir
->opts
.nFreeze
) {
1125 for(d
=0; d
<DIM
; d
++)
1126 if (pulldims
[d
] && ir
->opts
.nFreeze
[ggrpnr(groups
,egcFREEZE
,ii
)][d
])
1129 if (ir
->efep
== efepNO
) {
1132 m
= (1 - lambda
)*atom
->m
+ lambda
*atom
->mB
;
1134 if (pg
->nweight
> 0) {
1139 if (EI_ENERGY_MINIMIZATION(ir
->eI
)) {
1140 /* Move the mass to the weight */
1144 } else if (ir
->eI
== eiBD
) {
1146 mbd
= ir
->bd_fric
*ir
->delta_t
;
1148 if (groups
->grpnr
[egcTC
] == NULL
) {
1149 mbd
= ir
->delta_t
/ir
->opts
.tau_t
[0];
1151 mbd
= ir
->delta_t
/ir
->opts
.tau_t
[groups
->grpnr
[egcTC
][ii
]];
1163 gmx_mtop_atomlookup_destroy(alook
);
1166 gmx_fatal(FARGS
,"The total%s mass of pull group %d is zero",
1167 pg
->weight
? " weighted" : "",g
);
1171 "Pull group %d: %5d atoms, mass %9.3f",g
,pg
->nat
,tmass
);
1172 if (pg
->weight
|| EI_ENERGY_MINIMIZATION(ir
->eI
) || ir
->eI
== eiBD
) {
1173 fprintf(fplog
,", weighted mass %9.3f",wmass
*wmass
/wwmass
);
1175 if (pg
->epgrppbc
== epgrppbcCOS
) {
1176 fprintf(fplog
,", cosine weighting will be used");
1178 fprintf(fplog
,"\n");
1182 /* A value > 0 signals not frozen, it is updated later */
1186 for(d
=0; d
<DIM
; d
++)
1187 ndim
+= pulldims
[d
]*pg
->nat
;
1188 if (fplog
&& nfrozen
> 0 && nfrozen
< ndim
) {
1190 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1191 " that are subject to pulling are frozen.\n"
1192 " For pulling the whole group will be frozen.\n\n",
1200 void init_pull(FILE *fplog
,t_inputrec
*ir
,int nfile
,const t_filenm fnm
[],
1201 gmx_mtop_t
*mtop
,t_commrec
*cr
,const output_env_t oenv
, real lambda
,
1202 gmx_bool bOutFile
, unsigned long Flags
)
1206 int g
,start
=0,end
=0,m
;
1211 pull
->ePBC
= ir
->ePBC
;
1214 case epbcNONE
: pull
->npbcdim
= 0; break;
1215 case epbcXY
: pull
->npbcdim
= 2; break;
1216 default: pull
->npbcdim
= 3; break;
1221 fprintf(fplog
,"\nWill apply %s COM pulling in geometry '%s'\n",
1222 EPULLTYPE(ir
->ePull
),EPULLGEOM(pull
->eGeom
));
1223 if (pull
->grp
[0].nat
> 0)
1225 fprintf(fplog
,"between a reference group and %d group%s\n",
1226 pull
->ngrp
,pull
->ngrp
==1 ? "" : "s");
1230 fprintf(fplog
,"with an absolute reference on %d group%s\n",
1231 pull
->ngrp
,pull
->ngrp
==1 ? "" : "s");
1234 for(g
=0; g
<pull
->ngrp
+1; g
++)
1236 if (pull
->grp
[g
].nat
> 1 &&
1237 pull
->grp
[g
].pbcatom
< 0)
1239 /* We are using cosine weighting */
1240 fprintf(fplog
,"Cosine weighting is used for group %d\n",g
);
1246 please_cite(fplog
,"Engin2010");
1250 /* We always add the virial contribution,
1251 * except for geometry = direction_periodic where this is impossible.
1253 pull
->bVirial
= (pull
->eGeom
!= epullgDIRPBC
);
1254 if (getenv("GMX_NO_PULLVIR") != NULL
)
1258 fprintf(fplog
,"Found env. var., will not add the virial contribution of the COM pull forces\n");
1260 pull
->bVirial
= FALSE
;
1263 if (cr
&& PARTDECOMP(cr
))
1265 pd_at_range(cr
,&start
,&end
);
1269 pull
->dbuf_cyl
=NULL
;
1270 pull
->bRefAt
= FALSE
;
1272 for(g
=0; g
<pull
->ngrp
+1; g
++)
1274 pgrp
= &pull
->grp
[g
];
1275 pgrp
->epgrppbc
= epgrppbcNONE
;
1278 /* Determine if we need to take PBC into account for calculating
1279 * the COM's of the pull groups.
1281 for(m
=0; m
<pull
->npbcdim
; m
++)
1283 if (pull
->dim
[m
] && pgrp
->nat
> 1)
1285 if (pgrp
->pbcatom
>= 0)
1287 pgrp
->epgrppbc
= epgrppbcREFAT
;
1288 pull
->bRefAt
= TRUE
;
1294 gmx_fatal(FARGS
,"Pull groups can not have relative weights and cosine weighting at same time");
1296 pgrp
->epgrppbc
= epgrppbcCOS
;
1297 if (pull
->cosdim
>= 0 && pull
->cosdim
!= m
)
1299 gmx_fatal(FARGS
,"Can only use cosine weighting with pulling in one dimension (use mdp option pull_dim)");
1305 /* Set the indices */
1306 init_pull_group_index(fplog
,cr
,start
,end
,g
,pgrp
,pull
->dim
,mtop
,ir
,lambda
);
1307 if (PULL_CYL(pull
) && pgrp
->invtm
== 0)
1309 gmx_fatal(FARGS
,"Can not have frozen atoms in a cylinder pull group");
1314 /* Absolute reference, set the inverse mass to zero */
1320 /* if we use dynamic reference groups, do some initialising for them */
1323 if (pull
->grp
[0].nat
== 0)
1325 gmx_fatal(FARGS
, "Dynamic reference groups are not supported when using absolute reference!\n");
1327 snew(pull
->dyna
,pull
->ngrp
+1);
1330 /* Only do I/O when we are doing dynamics and if we are the MASTER */
1335 if (pull
->nstxout
> 0)
1337 pull
->out_x
= open_pull_out(opt2fn("-px",nfile
,fnm
),pull
,oenv
,TRUE
,Flags
);
1339 if (pull
->nstfout
> 0)
1341 pull
->out_f
= open_pull_out(opt2fn("-pf",nfile
,fnm
),pull
,oenv
,
1347 void finish_pull(FILE *fplog
,t_pull
*pull
)
1351 gmx_fio_fclose(pull
->out_x
);
1355 gmx_fio_fclose(pull
->out_f
);