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,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source 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.
44 #include "gromacs/domdec/localatomsetmanager.h"
45 #include "gromacs/fileio/readinp.h"
46 #include "gromacs/fileio/warninp.h"
47 #include "gromacs/gmxpreprocess/readir.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/mdlib/mdatoms.h"
50 #include "gromacs/mdtypes/inputrec.h"
51 #include "gromacs/mdtypes/md_enums.h"
52 #include "gromacs/mdtypes/pull_params.h"
53 #include "gromacs/pbcutil/pbc.h"
54 #include "gromacs/pulling/pull.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/cstringutil.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/futil.h"
59 #include "gromacs/utility/smalloc.h"
62 static void string2dvec(const char buf
[], dvec nums
)
66 if (sscanf(buf
, "%lf%lf%lf%lf", &nums
[0], &nums
[1], &nums
[2], &dum
) != 3)
68 gmx_fatal(FARGS
, "Expected three numbers at input line %s", buf
);
72 static void init_pull_group(t_pull_group
* pg
, const char* wbuf
)
78 while (sscanf(wbuf
, "%lf %n", &d
, &n
) == 1)
80 if (pg
->nweight
% 100 == 0)
82 srenew(pg
->weight
, pg
->nweight
+ 100);
84 pg
->weight
[pg
->nweight
++] = d
;
89 static void process_pull_dim(char* dim_buf
, ivec dim
, const t_pull_coord
* pcrd
)
92 char *ptr
, pulldim1
[STRLEN
];
96 for (d
= 0; d
< DIM
; d
++)
98 if (sscanf(ptr
, "%s%n", pulldim1
, &nchar
) != 1)
100 gmx_fatal(FARGS
, "Less than 3 pull dimensions given in pull_dim: '%s'", dim_buf
);
103 if (gmx::equalCaseInsensitive(pulldim1
, "N", 1))
107 else if (gmx::equalCaseInsensitive(pulldim1
, "Y", 1))
114 gmx_fatal(FARGS
, "Please use Y(ES) or N(O) for pull_dim only (not %s)", pulldim1
);
120 gmx_fatal(FARGS
, "All entries in pull dim are N");
122 if ((pcrd
->eGeom
== epullgDIHEDRAL
) && (ndim
< 3))
124 gmx_fatal(FARGS
, "Pull geometry dihedral is only useful with pull-dim = Y Y Y");
126 if ((pcrd
->eGeom
== epullgANGLE
|| pcrd
->eGeom
== epullgANGLEAXIS
) && (ndim
< 2))
129 "Pull geometry %s is only useful with pull-dim = Y for at least 2 dimensions",
130 EPULLGEOM(pcrd
->eGeom
));
134 static void init_pull_coord(t_pull_coord
* pcrd
,
135 int coord_index_for_output
,
137 const char* origin_buf
,
145 if (pcrd
->eType
== epullCONSTRAINT
146 && (pcrd
->eGeom
== epullgCYL
|| pcrd
->eGeom
== epullgDIRRELATIVE
|| pcrd
->eGeom
== epullgANGLE
147 || pcrd
->eGeom
== epullgANGLEAXIS
|| pcrd
->eGeom
== epullgDIHEDRAL
))
150 "Pulling of type %s can not be combined with geometry %s. Consider using pull "
152 epull_names
[pcrd
->eType
], epullg_names
[pcrd
->eGeom
], epull_names
[epullUMBRELLA
]);
155 if (pcrd
->eType
== epullEXTERNAL
)
157 if (pcrd
->externalPotentialProvider
[0] == '\0')
160 "The use of pull type '%s' for pull coordinate %d requires that the name of "
161 "the module providing the potential external is set with the option %s%d%s",
162 epull_names
[pcrd
->eType
], coord_index_for_output
, "pull-coord",
163 coord_index_for_output
, "-potential-provider");
164 warning_error(wi
, buf
);
170 "The use of pull type '%s' for pull coordinate %d requires that the pull rate "
172 epull_names
[pcrd
->eType
], coord_index_for_output
);
173 warning_error(wi
, buf
);
176 if (pcrd
->eGeom
== epullgCYL
)
178 /* Warn the user of a PBC restriction, caused by the fact that
179 * there is no reference value with an external pull potential.
182 "With pull type '%s' and geometry '%s', the distance component along the "
183 "cylinder axis between atoms in the cylinder group and the COM of the pull "
184 "group should be smaller than half the box length",
185 epull_names
[pcrd
->eType
], epullg_names
[pcrd
->eGeom
]);
186 warning_note(wi
, buf
);
190 process_pull_dim(dim_buf
, pcrd
->dim
, pcrd
);
192 string2dvec(origin_buf
, origin
);
193 if (pcrd
->group
[0] != 0 && dnorm(origin
) > 0)
195 gmx_fatal(FARGS
, "The pull origin can only be set with an absolute reference");
198 /* Check the given initial reference value and warn for dangerous values */
199 if (pcrd
->eGeom
== epullgDIST
)
201 if (pcrd
->bStart
&& pcrd
->init
< 0)
204 "The initial reference distance set by pull-coord-init is set to a negative "
205 "value (%g) with geometry %s while distances need to be non-negative. "
206 "This may work, since you have set pull-coord-start to 'yes' which modifies "
207 "this value, but only for certain starting distances. "
208 "If this is a mistake you may want to use geometry %s instead.",
209 pcrd
->init
, EPULLGEOM(pcrd
->eGeom
), EPULLGEOM(epullgDIR
));
213 else if (pcrd
->eGeom
== epullgANGLE
|| pcrd
->eGeom
== epullgANGLEAXIS
)
215 if (pcrd
->bStart
&& (pcrd
->init
< 0 || pcrd
->init
> 180))
217 /* This value of pcrd->init may be ok depending on pcrd->bStart which modifies pcrd->init later on */
219 "The initial reference angle set by pull-coord-init (%g) is outside of the "
220 "allowed range [0, 180] degrees for geometry (%s). "
221 "This may work, since you have set pull-coord-start to 'yes' which modifies "
222 "this value, but only for certain starting angles.",
223 pcrd
->init
, EPULLGEOM(pcrd
->eGeom
));
227 else if (pcrd
->eGeom
== epullgDIHEDRAL
)
229 if (pcrd
->bStart
&& (pcrd
->init
< -180 || pcrd
->init
> 180))
232 "The initial reference angle set by pull-coord-init (%g) is outside of the "
233 "allowed range [-180, 180] degrees for geometry (%s). "
234 "This may work, since you have set pull-coord-start to 'yes' which modifies "
235 "this value, but only for certain starting angles.",
236 pcrd
->init
, EPULLGEOM(pcrd
->eGeom
));
241 /* Check and set the pull vector */
243 string2dvec(vec_buf
, vec
);
245 if (pcrd
->eGeom
== epullgDIR
|| pcrd
->eGeom
== epullgCYL
|| pcrd
->eGeom
== epullgDIRPBC
246 || pcrd
->eGeom
== epullgANGLEAXIS
)
248 if (dnorm2(vec
) == 0)
250 gmx_fatal(FARGS
, "With pull geometry %s the pull vector can not be 0,0,0",
251 epullg_names
[pcrd
->eGeom
]);
253 for (int d
= 0; d
< DIM
; d
++)
255 if (vec
[d
] != 0 && pcrd
->dim
[d
] == 0)
258 "pull-coord-vec has non-zero %c-component while pull_dim for the "
259 "%c-dimension is set to N",
264 /* Normalize the direction vector */
265 dsvmul(1 / dnorm(vec
), vec
, vec
);
267 else /* This case is for are all the geometries where the pull vector is not used */
272 "A pull vector is given (%g %g %g) but will not be used with geometry %s. If "
273 "you really want to use this "
274 "vector, consider using geometry %s instead.",
275 vec
[0], vec
[1], vec
[2], EPULLGEOM(pcrd
->eGeom
),
276 pcrd
->eGeom
== epullgANGLE
? EPULLGEOM(epullgANGLEAXIS
) : EPULLGEOM(epullgDIR
));
280 for (m
= 0; m
< DIM
; m
++)
282 pcrd
->origin
[m
] = origin
[m
];
283 pcrd
->vec
[m
] = vec
[m
];
287 char** read_pullparams(std::vector
<t_inpfile
>* inp
, pull_params_t
* pull
, warninp_t wi
)
292 char provider
[STRLEN
], groups
[STRLEN
], dim_buf
[STRLEN
];
293 char wbuf
[STRLEN
], origin_buf
[STRLEN
], vec_buf
[STRLEN
];
298 /* read pull parameters */
299 printStringNoNewline(inp
, "Cylinder radius for dynamic reaction force groups (nm)");
300 pull
->cylinder_r
= get_ereal(inp
, "pull-cylinder-r", 1.5, wi
);
301 pull
->constr_tol
= get_ereal(inp
, "pull-constr-tol", 1E-6, wi
);
302 pull
->bPrintCOM
= (get_eeenum(inp
, "pull-print-com", yesno_names
, wi
) != 0);
303 pull
->bPrintRefValue
= (get_eeenum(inp
, "pull-print-ref-value", yesno_names
, wi
) != 0);
304 pull
->bPrintComp
= (get_eeenum(inp
, "pull-print-components", yesno_names
, wi
) != 0);
305 pull
->nstxout
= get_eint(inp
, "pull-nstxout", 50, wi
);
306 pull
->nstfout
= get_eint(inp
, "pull-nstfout", 50, wi
);
307 pull
->bSetPbcRefToPrevStepCOM
= (get_eeenum(inp
, "pull-pbc-ref-prev-step-com", yesno_names
, wi
) != 0);
308 pull
->bXOutAverage
= (get_eeenum(inp
, "pull-xout-average", yesno_names
, wi
) != 0);
309 pull
->bFOutAverage
= (get_eeenum(inp
, "pull-fout-average", yesno_names
, wi
) != 0);
310 printStringNoNewline(inp
, "Number of pull groups");
311 pull
->ngroup
= get_eint(inp
, "pull-ngroups", 1, wi
);
312 printStringNoNewline(inp
, "Number of pull coordinates");
313 pull
->ncoord
= get_eint(inp
, "pull-ncoords", 1, wi
);
315 if (pull
->ngroup
< 1)
317 gmx_fatal(FARGS
, "pull-ngroups should be >= 1");
319 /* We always add an absolute reference group (index 0), even if not used */
322 if (pull
->ncoord
< 1)
324 gmx_fatal(FARGS
, "pull-ncoords should be >= 1");
327 snew(pull
->group
, pull
->ngroup
);
329 snew(pull
->coord
, pull
->ncoord
);
331 /* pull group options */
332 printStringNoNewline(inp
, "Group and coordinate parameters");
334 /* Read the pull groups */
335 snew(grpbuf
, pull
->ngroup
);
336 /* Group 0 is the absolute reference, we don't read anything for 0 */
337 for (int groupNum
= 1; groupNum
< pull
->ngroup
; groupNum
++)
339 pgrp
= &pull
->group
[groupNum
];
340 snew(grpbuf
[groupNum
], STRLEN
);
341 sprintf(buf
, "pull-group%d-name", groupNum
);
342 setStringEntry(inp
, buf
, grpbuf
[groupNum
], "");
343 sprintf(buf
, "pull-group%d-weights", groupNum
);
344 setStringEntry(inp
, buf
, wbuf
, "");
345 sprintf(buf
, "pull-group%d-pbcatom", groupNum
);
346 pgrp
->pbcatom
= get_eint(inp
, buf
, 0, wi
);
348 /* Initialize the pull group */
349 init_pull_group(pgrp
, wbuf
);
352 /* Read the pull coordinates */
353 for (int coordNum
= 1; coordNum
< pull
->ncoord
+ 1; coordNum
++)
355 pcrd
= &pull
->coord
[coordNum
- 1];
356 sprintf(buf
, "pull-coord%d-type", coordNum
);
357 pcrd
->eType
= get_eeenum(inp
, buf
, epull_names
, wi
);
358 sprintf(buf
, "pull-coord%d-potential-provider", coordNum
);
359 setStringEntry(inp
, buf
, provider
, "");
360 pcrd
->externalPotentialProvider
= gmx_strdup(provider
);
361 sprintf(buf
, "pull-coord%d-geometry", coordNum
);
362 pcrd
->eGeom
= get_eeenum(inp
, buf
, epullg_names
, wi
);
363 sprintf(buf
, "pull-coord%d-groups", coordNum
);
364 setStringEntry(inp
, buf
, groups
, "");
368 case epullgDIHEDRAL
: pcrd
->ngroup
= 6; break;
369 case epullgDIRRELATIVE
:
370 case epullgANGLE
: pcrd
->ngroup
= 4; break;
371 default: pcrd
->ngroup
= 2; break;
374 nscan
= sscanf(groups
, "%d %d %d %d %d %d %d", &pcrd
->group
[0], &pcrd
->group
[1],
375 &pcrd
->group
[2], &pcrd
->group
[3], &pcrd
->group
[4], &pcrd
->group
[5], &idum
);
376 if (nscan
!= pcrd
->ngroup
)
379 gmx::formatString("%s should contain %d pull group indices with geometry %s",
380 buf
, pcrd
->ngroup
, epullg_names
[pcrd
->eGeom
]);
381 set_warning_line(wi
, nullptr, -1);
382 warning_error(wi
, message
);
384 for (int g
= 0; g
< pcrd
->ngroup
; g
++)
386 if (pcrd
->group
[g
] < 0 || pcrd
->group
[g
] >= pull
->ngroup
)
388 /* Quit with a fatal error to avoid invalid memory access */
390 "%s contains an invalid pull group %d, you should have %d <= group <= %d",
391 buf
, pcrd
->group
[g
], 0, pull
->ngroup
- 1);
395 sprintf(buf
, "pull-coord%d-dim", coordNum
);
396 setStringEntry(inp
, buf
, dim_buf
, "Y Y Y");
397 sprintf(buf
, "pull-coord%d-origin", coordNum
);
398 setStringEntry(inp
, buf
, origin_buf
, "0.0 0.0 0.0");
399 sprintf(buf
, "pull-coord%d-vec", coordNum
);
400 setStringEntry(inp
, buf
, vec_buf
, "0.0 0.0 0.0");
401 sprintf(buf
, "pull-coord%d-start", coordNum
);
402 pcrd
->bStart
= (get_eeenum(inp
, buf
, yesno_names
, wi
) != 0);
403 sprintf(buf
, "pull-coord%d-init", coordNum
);
404 pcrd
->init
= get_ereal(inp
, buf
, 0.0, wi
);
405 sprintf(buf
, "pull-coord%d-rate", coordNum
);
406 pcrd
->rate
= get_ereal(inp
, buf
, 0.0, wi
);
407 sprintf(buf
, "pull-coord%d-k", coordNum
);
408 pcrd
->k
= get_ereal(inp
, buf
, 0.0, wi
);
409 sprintf(buf
, "pull-coord%d-kB", coordNum
);
410 pcrd
->kB
= get_ereal(inp
, buf
, pcrd
->k
, wi
);
412 /* Initialize the pull coordinate */
413 init_pull_coord(pcrd
, coordNum
, dim_buf
, origin_buf
, vec_buf
, wi
);
419 void make_pull_groups(pull_params_t
* pull
, char** pgnames
, const t_blocka
* grps
, char** gnames
)
424 /* Absolute reference group (might not be used) is special */
425 pgrp
= &pull
->group
[0];
428 pgrp
->pbcatom_input
= -1;
430 for (g
= 1; g
< pull
->ngroup
; g
++)
432 pgrp
= &pull
->group
[g
];
433 pgrp
->pbcatom_input
= pgrp
->pbcatom
;
435 if (strcmp(pgnames
[g
], "") == 0)
437 gmx_fatal(FARGS
, "Pull option pull_group%d required by grompp has not been set.", g
);
440 ig
= search_string(pgnames
[g
], grps
->nr
, gnames
);
441 pgrp
->nat
= grps
->index
[ig
+ 1] - grps
->index
[ig
];
443 fprintf(stderr
, "Pull group %d '%s' has %d atoms\n", g
, pgnames
[g
], pgrp
->nat
);
447 gmx_fatal(FARGS
, "Pull group %d '%s' is empty", g
, pgnames
[g
]);
450 snew(pgrp
->ind
, pgrp
->nat
);
451 for (i
= 0; i
< pgrp
->nat
; i
++)
453 pgrp
->ind
[i
] = grps
->a
[grps
->index
[ig
] + i
];
456 if (pgrp
->nweight
> 0 && pgrp
->nweight
!= pgrp
->nat
)
459 "Number of weights (%d) for pull group %d '%s' does not match the number of "
461 pgrp
->nweight
, g
, pgnames
[g
], pgrp
->nat
);
466 /* No pbc is required for this group */
471 if (pgrp
->pbcatom
> 0)
475 else if (pgrp
->pbcatom
== 0)
477 pgrp
->pbcatom
= pgrp
->ind
[(pgrp
->nat
- 1) / 2];
481 /* Use cosine weighting */
488 void make_pull_coords(pull_params_t
* pull
)
493 for (c
= 0; c
< pull
->ncoord
; c
++)
495 pcrd
= &pull
->coord
[c
];
497 if (pcrd
->group
[0] < 0 || pcrd
->group
[0] >= pull
->ngroup
|| pcrd
->group
[1] < 0
498 || pcrd
->group
[1] >= pull
->ngroup
)
501 "Pull group index in pull-coord%d-groups out of range, should be between %d "
503 c
+ 1, 0, pull
->ngroup
+ 1);
506 if (pcrd
->group
[0] == pcrd
->group
[1])
508 gmx_fatal(FARGS
, "Identical pull group indices in pull-coord%d-groups", c
+ 1);
511 if (pcrd
->eGeom
== epullgCYL
)
513 if (pull
->group
[pcrd
->group
[0]].nweight
> 0)
517 "Weights are not supported for the reference group with cylinder pulling");
523 pull_t
* set_pull_init(t_inputrec
* ir
, const gmx_mtop_t
* mtop
, rvec
* x
, matrix box
, real lambda
, warninp_t wi
)
532 gmx::LocalAtomSetManager atomSets
;
533 pull_work
= init_pull(nullptr, pull
, ir
, mtop
, nullptr, &atomSets
, lambda
);
534 auto mdAtoms
= gmx::makeMDAtoms(nullptr, *mtop
, *ir
, false);
535 auto md
= mdAtoms
->mdatoms();
536 atoms2md(mtop
, ir
, -1, nullptr, mtop
->natoms
, mdAtoms
.get());
539 update_mdatoms(md
, lambda
);
542 set_pbc(&pbc
, ir
->ePBC
, box
);
544 t_start
= ir
->init_t
+ ir
->init_step
* ir
->delta_t
;
546 if (pull
->bSetPbcRefToPrevStepCOM
)
548 initPullComFromPrevStep(nullptr, pull_work
, md
, &pbc
, x
);
550 pull_calc_coms(nullptr, pull_work
, md
, &pbc
, t_start
, x
, nullptr);
552 for (int g
= 0; g
< pull
->ngroup
; g
++)
554 bool groupObeysPbc
= pullCheckPbcWithinGroup(
555 *pull_work
, gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec
*>(x
), mtop
->natoms
),
556 pbc
, g
, c_pullGroupSmallGroupThreshold
);
560 if (pull
->group
[g
].pbcatom_input
== 0)
563 "When the maximum distance from a pull group reference atom to other atoms "
565 "group is larger than %g times half the box size a centrally placed "
566 "atom should be chosen as pbcatom. Pull group %d is larger than that and "
568 "a specific atom selected as reference atom.",
569 c_pullGroupSmallGroupThreshold
, g
);
570 warning_error(wi
, buf
);
572 else if (!pull
->bSetPbcRefToPrevStepCOM
)
575 "The maximum distance from the chosen PBC atom (%d) of pull group %d to "
577 "atoms in the group is larger than %g times half the box size. "
578 "Set the pull-pbc-ref-prev-step-com option to yes.",
579 pull
->group
[g
].pbcatom
+ 1, g
, c_pullGroupSmallGroupThreshold
);
580 warning_error(wi
, buf
);
585 groupObeysPbc
= pullCheckPbcWithinGroup(
586 *pull_work
, gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec
*>(x
), mtop
->natoms
),
587 pbc
, g
, c_pullGroupPbcMargin
);
592 "Pull group %d has atoms at a distance larger than %g times half the box "
593 "size from the PBC atom (%d). "
594 "If atoms are or will more beyond half the box size from the PBC atom, the "
595 "COM will be ill defined.",
596 g
, c_pullGroupPbcMargin
, pull
->group
[g
].pbcatom
+ 1);
597 set_warning_line(wi
, nullptr, -1);
603 fprintf(stderr
, "Pull group natoms pbc atom distance at start reference at t=0\n");
604 for (c
= 0; c
< pull
->ncoord
; c
++)
607 t_pull_group
*pgrp0
, *pgrp1
;
611 pcrd
= &pull
->coord
[c
];
613 pgrp0
= &pull
->group
[pcrd
->group
[0]];
614 pgrp1
= &pull
->group
[pcrd
->group
[1]];
615 fprintf(stderr
, "%8d %8d %8d\n", pcrd
->group
[0], pgrp0
->nat
, pgrp0
->pbcatom
+ 1);
616 fprintf(stderr
, "%8d %8d %8d ", pcrd
->group
[1], pgrp1
->nat
, pgrp1
->pbcatom
+ 1);
624 value
= get_pull_coord_value(pull_work
, c
, &pbc
);
626 value
*= pull_conversion_factor_internal2userinput(pcrd
);
627 fprintf(stderr
, " %10.3f %s", value
, pull_coordinate_units(pcrd
));
631 pcrd
->init
= value
+ init
;
634 if (pcrd
->eGeom
== epullgDIST
)
639 "The initial pull distance (%g) needs to be non-negative with geometry "
640 "%s. If you want a signed distance, use geometry %s instead.",
641 pcrd
->init
, EPULLGEOM(pcrd
->eGeom
), EPULLGEOM(epullgDIR
));
644 /* TODO: With a positive init but a negative rate things could still
645 * go wrong, but it might be fine if you don't pull too far.
646 * We should give a warning or note when there is only one pull dim
647 * active, since that is usually the problematic case when you should
648 * be using direction. We will do this later, since an already planned
649 * generalization of the pull code makes pull dim available here.
652 else if (pcrd
->eGeom
== epullgANGLE
|| pcrd
->eGeom
== epullgANGLEAXIS
)
654 if (pcrd
->init
< 0 || pcrd
->init
> 180)
657 "The initial pull reference angle (%g) is outside of the allowed range "
662 else if (pcrd
->eGeom
== epullgDIHEDRAL
)
664 if (pcrd
->init
< -180 || pcrd
->init
> 180)
667 "The initial pull reference angle (%g) is outside of the allowed range "
668 "[-180, 180] degrees.",
674 fprintf(stderr
, " %10.3f %s\n", pcrd
->init
, pull_coordinate_units(pcrd
));