3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
49 #include "gmx_fatal.h"
60 static char pulldim
[STRLEN
];
62 static void string2dvec(char buf
[], dvec nums
)
64 if (sscanf(buf
,"%lf%lf%lf",&nums
[0],&nums
[1],&nums
[2]) != 3)
65 gmx_fatal(FARGS
,"Expected three numbers at input line %s",buf
);
68 static void init_pullgrp(t_pullgrp
*pg
,char *wbuf
,
69 bool bRef
,int eGeom
,char *s_vec
)
76 while (sscanf(wbuf
,"%lf %n",&d
,&n
) == 1) {
77 if (pg
->nweight
% 100 == 0) {
78 srenew(pg
->weight
,pg
->nweight
+100);
80 pg
->weight
[pg
->nweight
++] = d
;
84 if (eGeom
== epullgDIST
) {
87 string2dvec(s_vec
,vec
);
88 if (eGeom
== epullgDIR
|| eGeom
== epullgCYL
||
89 (eGeom
== epullgPOS
&& dnorm(vec
) != 0))
90 /* Normalize the direction vector */
91 dsvmul(1/dnorm(vec
),vec
,vec
);
98 char **read_pullparams(int *ninp_p
,t_inpfile
**inp_p
,
99 t_pull
*pull
,bool *bStart
,int *nerror_p
)
101 int ninp
,nerror
,i
,nchar
,ndim
,nscan
,m
;
105 char dummy
[STRLEN
],buf
[STRLEN
],init
[STRLEN
];
106 const char *init_def1
="0.0",*init_def3
="0.0 0.0 0.0";
107 char wbuf
[STRLEN
],VecTemp
[STRLEN
];
116 /* read pull parameters */
117 CTYPE("Pull geometry: distance, direction, cylinder or position");
118 EETYPE("pull_geometry", pull
->eGeom
, epullg_names
, &nerror
, TRUE
);
119 CTYPE("Select components for the pull vector. default: Y Y Y");
120 STYPE("pull_dim", pulldim
, "Y Y Y");
121 CTYPE("Cylinder radius for dynamic reaction force groups (nm)");
122 RTYPE("pull_r1", pull
->cyl_r1
, 1.0);
123 CTYPE("Switch from r1 to r0 in case of dynamic reaction force");
124 RTYPE("pull_r0", pull
->cyl_r0
, 1.5);
125 RTYPE("pull_constr_tol", pull
->constr_tol
, 1E-6);
126 EETYPE("pull_start", *bStart
, yesno_names
, &nerror
, TRUE
);
127 ITYPE("pull_nstxout", pull
->nstxout
, 10);
128 ITYPE("pull_nstfout", pull
->nstfout
, 1);
129 CTYPE("Number of pull groups");
130 ITYPE("pull_ngroups", pull
->ngrp
,1);
132 if (pull
->cyl_r1
> pull
->cyl_r0
) {
133 fprintf(stderr
,"ERROR: pull_r1 > pull_r0\n");
137 if (pull
->ngrp
< 1) {
138 gmx_fatal(FARGS
,"pull_ngroups should be >= 1");
141 snew(pull
->grp
,pull
->ngrp
+1);
143 if (pull
->eGeom
== epullgPOS
) {
149 /* pull group options */
150 CTYPE("Group name, weight (default all 1), vector, init, rate (nm/ps), kJ/(mol*nm^2)");
151 /* Read the pull groups */
152 snew(grpbuf
,pull
->ngrp
+1);
153 for(i
=0; i
<pull
->ngrp
+1; i
++) {
154 pgrp
= &pull
->grp
[i
];
155 snew(grpbuf
[i
],STRLEN
);
156 sprintf(buf
,"pull_group%d",i
);
157 STYPE(buf
, grpbuf
[i
], "");
158 sprintf(buf
,"pull_weights%d",i
);
159 STYPE(buf
, wbuf
, "");
160 sprintf(buf
,"pull_pbcatom%d",i
);
161 ITYPE(buf
, pgrp
->pbcatom
, 0);
163 sprintf(buf
,"pull_vec%d",i
);
164 STYPE(buf
, VecTemp
, "0.0 0.0 0.0");
165 sprintf(buf
,"pull_init%d",i
);
166 STYPE(buf
, init
, ndim
==1 ? init_def1
: init_def3
);
167 nscan
= sscanf(init
,"%lf %lf %lf",&vec
[0],&vec
[1],&vec
[2]);
169 fprintf(stderr
,"ERROR: %s should have %d components\n",buf
,ndim
);
172 for(m
=0; m
<DIM
; m
++) {
173 pgrp
->init
[m
] = (m
<ndim
? vec
[m
] : 0.0);
175 sprintf(buf
,"pull_rate%d",i
);
176 RTYPE(buf
, pgrp
->rate
, 0.0);
177 sprintf(buf
,"pull_k%d",i
);
178 RTYPE(buf
, pgrp
->k
, 0.0);
179 sprintf(buf
,"pull_kB%d",i
);
180 RTYPE(buf
, pgrp
->kB
, pgrp
->k
);
183 /* Initialize the pull group */
184 init_pullgrp(pgrp
,wbuf
,i
==0,pull
->eGeom
,VecTemp
);
194 void make_pull_groups(t_pull
*pull
,char **pgnames
,t_blocka
*grps
,char **gnames
)
196 int d
,nchar
,g
,ig
=-1,i
;
197 char *ptr
,pulldim1
[STRLEN
];
202 for(d
=0; d
<DIM
; d
++) {
203 if (sscanf(ptr
,"%s%n",pulldim1
,&nchar
) != 1)
204 gmx_fatal(FARGS
,"Less than 3 pull dimensions given in pull_dim: '%s'",
207 if (strncasecmp(pulldim1
,"N",1) == 0) {
209 } else if (strncasecmp(pulldim1
,"Y",1) == 0) {
213 gmx_fatal(FARGS
,"Please use Y(ES) or N(O) for pull_dim only (not %s)",
219 gmx_fatal(FARGS
,"All entries in pull_dim are N");
221 for(g
=0; g
<pull
->ngrp
+1; g
++) {
222 pgrp
= &pull
->grp
[g
];
223 if (g
== 0 && strcmp(pgnames
[g
],"") == 0) {
226 ig
= search_string(pgnames
[g
],grps
->nr
,gnames
);
227 pgrp
->nat
= grps
->index
[ig
+1] - grps
->index
[ig
];
230 fprintf(stderr
,"Pull group %d '%s' has %d atoms\n",
231 g
,pgnames
[g
],pgrp
->nat
);
232 snew(pgrp
->ind
,pgrp
->nat
);
233 for(i
=0; i
<pgrp
->nat
; i
++)
234 pgrp
->ind
[i
] = grps
->a
[grps
->index
[ig
]+i
];
236 if (pull
->eGeom
== epullgCYL
&& g
== 0 && pgrp
->nweight
> 0)
237 gmx_fatal(FARGS
,"Weights are not supported for the reference group with cylinder pulling");
238 if (pgrp
->nweight
> 0 && pgrp
->nweight
!= pgrp
->nat
)
239 gmx_fatal(FARGS
,"Number of weights (%d) for pull group %d '%s' does not match the number of atoms (%d)",
240 pgrp
->nweight
,g
,pgnames
[g
],pgrp
->nat
);
242 if (pgrp
->nat
== 1) {
243 /* No pbc is required for this group */
246 if (pgrp
->pbcatom
> 0) {
248 } else if (pgrp
->pbcatom
== 0) {
249 pgrp
->pbcatom
= pgrp
->ind
[(pgrp
->nat
-1)/2];
251 /* Use cosine weighting */
256 if (g
> 0 && pull
->eGeom
!= epullgDIST
) {
257 for(d
=0; d
<DIM
; d
++) {
258 if (pgrp
->vec
[d
] != 0 && pull
->dim
[d
] == 0) {
259 gmx_fatal(FARGS
,"ERROR: pull_vec%d has non-zero %c-component while pull_dim in N\n",g
,'x'+d
);
264 if ((pull
->eGeom
== epullgDIR
|| pull
->eGeom
== epullgCYL
) &&
265 g
> 0 && norm2(pgrp
->vec
) == 0)
266 gmx_fatal(FARGS
,"pull_vec%d can not be zero with geometry %s",
267 g
,EPULLGEOM(pull
->eGeom
));
268 if ((pull
->eGeom
== epullgPOS
) && pgrp
->rate
!= 0 &&
269 g
> 0 && norm2(pgrp
->vec
) == 0)
270 gmx_fatal(FARGS
,"pull_vec%d can not be zero with geometry %s and non-zero rate",
271 g
,EPULLGEOM(pull
->eGeom
));
274 if (pull
->eGeom
== epullgCYL
)
275 gmx_fatal(FARGS
,"Absolute reference groups are not supported with geometry %s",EPULLGEOM(pull
->eGeom
));
277 gmx_fatal(FARGS
,"Pull group %d '%s' is empty",g
,pgnames
[g
]);
284 void set_pull_init(t_inputrec
*ir
,gmx_mtop_t
*mtop
,rvec
*x
,matrix box
,
285 const output_env_t oenv
,bool bStart
)
292 double t_start
,tinvrate
;
296 init_pull(NULL
,ir
,0,NULL
,mtop
,NULL
,oenv
,FALSE
,0);
297 md
= init_mdatoms(NULL
,mtop
,ir
->efep
);
298 atoms2md(mtop
,ir
,0,NULL
,0,mtop
->natoms
,md
);
300 update_mdatoms(md
,ir
->init_lambda
);
303 if (pull
->eGeom
== epullgPOS
)
308 set_pbc(&pbc
,ir
->ePBC
,box
);
310 t_start
= ir
->init_t
+ ir
->init_step
*ir
->delta_t
;
312 pull_calc_coms(NULL
,pull
,md
,&pbc
,t_start
,x
,NULL
);
314 fprintf(stderr
,"Pull group natoms pbc atom distance at start reference at t=0\n");
315 for(g
=0; g
<pull
->ngrp
+1; g
++) {
316 pgrp
= &pull
->grp
[g
];
317 fprintf(stderr
,"%8d %8d %8d ",g
,pgrp
->nat
,pgrp
->pbcatom
+1);
318 copy_rvec(pgrp
->init
,init
);
319 clear_rvec(pgrp
->init
);
324 tinvrate
= t_start
/pgrp
->rate
;
325 get_pullgrp_distance(pull
,&pbc
,g
,0,dr
,dev
);
326 for(m
=0; m
<DIM
; m
++) {
328 fprintf(stderr
," %6.3f",dev
[m
]);
333 for(m
=0; m
<DIM
; m
++) {
335 pgrp
->init
[m
] = init
[m
] + dev
[m
]
336 - tinvrate
*(pull
->eGeom
==epullgPOS
? pgrp
->vec
[m
] : 1);
338 pgrp
->init
[m
] = init
[m
];
340 fprintf(stderr
," %6.3f",pgrp
->init
[m
]);
345 fprintf(stderr
,"\n");