File: | gromacs/gmxana/gmx_trjconv.c |
Location: | line 1332, column 21 |
Description: | Function call argument is an uninitialized value |
1 | /* | |||
2 | * This file is part of the GROMACS molecular simulation package. | |||
3 | * | |||
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. | |||
5 | * Copyright (c) 2001-2004, The GROMACS development team. | |||
6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by | |||
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, | |||
8 | * and including many others, as listed in the AUTHORS file in the | |||
9 | * top-level source directory and at http://www.gromacs.org. | |||
10 | * | |||
11 | * GROMACS is free software; you can redistribute it and/or | |||
12 | * modify it under the terms of the GNU Lesser General Public License | |||
13 | * as published by the Free Software Foundation; either version 2.1 | |||
14 | * of the License, or (at your option) any later version. | |||
15 | * | |||
16 | * GROMACS is distributed in the hope that it will be useful, | |||
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |||
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |||
19 | * Lesser General Public License for more details. | |||
20 | * | |||
21 | * You should have received a copy of the GNU Lesser General Public | |||
22 | * License along with GROMACS; if not, see | |||
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, | |||
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. | |||
25 | * | |||
26 | * If you want to redistribute modifications to GROMACS, please | |||
27 | * consider that scientific software is very special. Version | |||
28 | * control is crucial - bugs must be traceable. We will be happy to | |||
29 | * consider code for inclusion in the official distribution, but | |||
30 | * derived work must not be called official GROMACS. Details are found | |||
31 | * in the README & COPYING files - if they are missing, get the | |||
32 | * official version at http://www.gromacs.org. | |||
33 | * | |||
34 | * To help us fund GROMACS development, we humbly ask that you cite | |||
35 | * the research papers on the package. Check out http://www.gromacs.org. | |||
36 | */ | |||
37 | #ifdef HAVE_CONFIG_H1 | |||
38 | #include <config.h> | |||
39 | #endif | |||
40 | ||||
41 | #include <math.h> | |||
42 | #include <stdlib.h> | |||
43 | #include <string.h> | |||
44 | ||||
45 | #include "copyrite.h" | |||
46 | #include "macros.h" | |||
47 | #include "gromacs/utility/smalloc.h" | |||
48 | #include "typedefs.h" | |||
49 | #include "gromacs/fileio/gmxfio.h" | |||
50 | #include "gromacs/fileio/tpxio.h" | |||
51 | #include "gromacs/fileio/trxio.h" | |||
52 | #include "gromacs/fileio/trnio.h" | |||
53 | #include "gromacs/fileio/tngio_for_tools.h" | |||
54 | #include "gromacs/commandline/pargs.h" | |||
55 | #include "gromacs/utility/futil.h" | |||
56 | #include "gromacs/fileio/pdbio.h" | |||
57 | #include "gromacs/fileio/confio.h" | |||
58 | #include "names.h" | |||
59 | #include "index.h" | |||
60 | #include "gromacs/math/vec.h" | |||
61 | #include "gromacs/fileio/xtcio.h" | |||
62 | #include "rmpbc.h" | |||
63 | #include "pbc.h" | |||
64 | #include "viewit.h" | |||
65 | #include "gromacs/fileio/xvgr.h" | |||
66 | #include "gmx_ana.h" | |||
67 | ||||
68 | #include "gromacs/math/do_fit.h" | |||
69 | #include "gromacs/utility/fatalerror.h" | |||
70 | ||||
71 | #ifdef HAVE_UNISTD_H | |||
72 | #include <unistd.h> | |||
73 | #endif | |||
74 | ||||
75 | enum { | |||
76 | euSel, euRect, euTric, euCompact, euNR | |||
77 | }; | |||
78 | ||||
79 | ||||
80 | static void calc_pbc_cluster(int ecenter, int nrefat, t_topology *top, int ePBC, | |||
81 | rvec x[], atom_id index[], matrix box) | |||
82 | { | |||
83 | int m, i, j, j0, j1, jj, ai, aj; | |||
84 | int imin, jmin; | |||
85 | real fac, min_dist2; | |||
86 | rvec dx, xtest, box_center; | |||
87 | int nmol, imol_center; | |||
88 | atom_id *molind; | |||
89 | gmx_bool *bMol, *bTmp; | |||
90 | rvec *m_com, *m_shift; | |||
91 | t_pbc pbc; | |||
92 | real *com_dist2; | |||
93 | int *cluster; | |||
94 | int *added; | |||
95 | int ncluster, nadded; | |||
96 | real tmp_r2; | |||
97 | ||||
98 | calc_box_center(ecenter, box, box_center); | |||
99 | ||||
100 | /* Initiate the pbc structure */ | |||
101 | memset(&pbc, 0, sizeof(pbc)); | |||
102 | set_pbc(&pbc, ePBC, box); | |||
103 | ||||
104 | /* Convert atom index to molecular */ | |||
105 | nmol = top->mols.nr; | |||
106 | molind = top->mols.index; | |||
107 | snew(bMol, nmol)(bMol) = save_calloc("bMol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 107, (nmol), sizeof(*(bMol))); | |||
108 | snew(m_com, nmol)(m_com) = save_calloc("m_com", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 108, (nmol), sizeof(*(m_com))); | |||
109 | snew(m_shift, nmol)(m_shift) = save_calloc("m_shift", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 109, (nmol), sizeof(*(m_shift))); | |||
110 | snew(cluster, nmol)(cluster) = save_calloc("cluster", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 110, (nmol), sizeof(*(cluster))); | |||
111 | snew(added, nmol)(added) = save_calloc("added", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 111, (nmol), sizeof(*(added))); | |||
112 | snew(bTmp, top->atoms.nr)(bTmp) = save_calloc("bTmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 112, (top->atoms.nr), sizeof(*(bTmp))); | |||
113 | ||||
114 | for (i = 0; (i < nrefat); i++) | |||
115 | { | |||
116 | /* Mark all molecules in the index */ | |||
117 | ai = index[i]; | |||
118 | bTmp[ai] = TRUE1; | |||
119 | /* Binary search assuming the molecules are sorted */ | |||
120 | j0 = 0; | |||
121 | j1 = nmol-1; | |||
122 | while (j0 < j1) | |||
123 | { | |||
124 | if (ai < molind[j0+1]) | |||
125 | { | |||
126 | j1 = j0; | |||
127 | } | |||
128 | else if (ai >= molind[j1]) | |||
129 | { | |||
130 | j0 = j1; | |||
131 | } | |||
132 | else | |||
133 | { | |||
134 | jj = (j0+j1)/2; | |||
135 | if (ai < molind[jj+1]) | |||
136 | { | |||
137 | j1 = jj; | |||
138 | } | |||
139 | else | |||
140 | { | |||
141 | j0 = jj; | |||
142 | } | |||
143 | } | |||
144 | } | |||
145 | bMol[j0] = TRUE1; | |||
146 | } | |||
147 | /* Double check whether all atoms in all molecules that are marked are part | |||
148 | * of the cluster. Simultaneously compute the center of geometry. | |||
149 | */ | |||
150 | min_dist2 = 10*sqr(trace(box)); | |||
151 | imol_center = -1; | |||
152 | ncluster = 0; | |||
153 | for (i = 0; i < nmol; i++) | |||
154 | { | |||
155 | for (j = molind[i]; j < molind[i+1]; j++) | |||
156 | { | |||
157 | if (bMol[i] && !bTmp[j]) | |||
158 | { | |||
159 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 159, "Molecule %d marked for clustering but not atom %d in it - check your index!", i+1, j+1); | |||
160 | } | |||
161 | else if (!bMol[i] && bTmp[j]) | |||
162 | { | |||
163 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 163, "Atom %d marked for clustering but not molecule %d - this is an internal error...", j+1, i+1); | |||
164 | } | |||
165 | else if (bMol[i]) | |||
166 | { | |||
167 | /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */ | |||
168 | if (j > molind[i]) | |||
169 | { | |||
170 | pbc_dx(&pbc, x[j], x[j-1], dx); | |||
171 | rvec_add(x[j-1], dx, x[j]); | |||
172 | } | |||
173 | /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */ | |||
174 | rvec_inc(m_com[i], x[j]); | |||
175 | } | |||
176 | } | |||
177 | if (bMol[i]) | |||
178 | { | |||
179 | /* Normalize center of geometry */ | |||
180 | fac = 1.0/(molind[i+1]-molind[i]); | |||
181 | for (m = 0; (m < DIM3); m++) | |||
182 | { | |||
183 | m_com[i][m] *= fac; | |||
184 | } | |||
185 | /* Determine which molecule is closest to the center of the box */ | |||
186 | pbc_dx(&pbc, box_center, m_com[i], dx); | |||
187 | tmp_r2 = iprod(dx, dx); | |||
188 | ||||
189 | if (tmp_r2 < min_dist2) | |||
190 | { | |||
191 | min_dist2 = tmp_r2; | |||
192 | imol_center = i; | |||
193 | } | |||
194 | cluster[ncluster++] = i; | |||
195 | } | |||
196 | } | |||
197 | sfree(bTmp)save_free("bTmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 197, (bTmp)); | |||
198 | ||||
199 | if (ncluster <= 0) | |||
200 | { | |||
201 | fprintf(stderrstderr, "No molecules selected in the cluster\n"); | |||
202 | return; | |||
203 | } | |||
204 | else if (imol_center == -1) | |||
205 | { | |||
206 | fprintf(stderrstderr, "No central molecules could be found\n"); | |||
207 | return; | |||
208 | } | |||
209 | ||||
210 | nadded = 0; | |||
211 | added[nadded++] = imol_center; | |||
212 | bMol[imol_center] = FALSE0; | |||
213 | ||||
214 | while (nadded < ncluster) | |||
215 | { | |||
216 | /* Find min distance between cluster molecules and those remaining to be added */ | |||
217 | min_dist2 = 10*sqr(trace(box)); | |||
218 | imin = -1; | |||
219 | jmin = -1; | |||
220 | /* Loop over added mols */ | |||
221 | for (i = 0; i < nadded; i++) | |||
222 | { | |||
223 | ai = added[i]; | |||
224 | /* Loop over all mols */ | |||
225 | for (j = 0; j < ncluster; j++) | |||
226 | { | |||
227 | aj = cluster[j]; | |||
228 | /* check those remaining to be added */ | |||
229 | if (bMol[aj]) | |||
230 | { | |||
231 | pbc_dx(&pbc, m_com[aj], m_com[ai], dx); | |||
232 | tmp_r2 = iprod(dx, dx); | |||
233 | if (tmp_r2 < min_dist2) | |||
234 | { | |||
235 | min_dist2 = tmp_r2; | |||
236 | imin = ai; | |||
237 | jmin = aj; | |||
238 | } | |||
239 | } | |||
240 | } | |||
241 | } | |||
242 | ||||
243 | /* Add the best molecule */ | |||
244 | added[nadded++] = jmin; | |||
245 | bMol[jmin] = FALSE0; | |||
246 | /* Calculate the shift from the ai molecule */ | |||
247 | pbc_dx(&pbc, m_com[jmin], m_com[imin], dx); | |||
248 | rvec_add(m_com[imin], dx, xtest); | |||
249 | rvec_sub(xtest, m_com[jmin], m_shift[jmin]); | |||
250 | rvec_inc(m_com[jmin], m_shift[jmin]); | |||
251 | ||||
252 | for (j = molind[jmin]; j < molind[jmin+1]; j++) | |||
253 | { | |||
254 | rvec_inc(x[j], m_shift[jmin]); | |||
255 | } | |||
256 | fprintf(stdoutstdout, "\rClustering iteration %d of %d...", nadded, ncluster); | |||
257 | fflush(stdoutstdout); | |||
258 | } | |||
259 | ||||
260 | sfree(added)save_free("added", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 260, (added)); | |||
261 | sfree(cluster)save_free("cluster", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 261, (cluster)); | |||
262 | sfree(bMol)save_free("bMol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 262, (bMol)); | |||
263 | sfree(m_com)save_free("m_com", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 263, (m_com)); | |||
264 | sfree(m_shift)save_free("m_shift", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 264, (m_shift)); | |||
265 | ||||
266 | fprintf(stdoutstdout, "\n"); | |||
267 | } | |||
268 | ||||
269 | static void put_molecule_com_in_box(int unitcell_enum, int ecenter, | |||
270 | t_block *mols, | |||
271 | int natoms, t_atom atom[], | |||
272 | int ePBC, matrix box, rvec x[]) | |||
273 | { | |||
274 | atom_id i, j; | |||
275 | int d; | |||
276 | rvec com, new_com, shift, dx, box_center; | |||
277 | real m; | |||
278 | double mtot; | |||
279 | t_pbc pbc; | |||
280 | ||||
281 | calc_box_center(ecenter, box, box_center); | |||
282 | set_pbc(&pbc, ePBC, box); | |||
283 | if (mols->nr <= 0) | |||
284 | { | |||
285 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 285, "There are no molecule descriptions. I need a .tpr file for this pbc option."); | |||
286 | } | |||
287 | for (i = 0; (i < mols->nr); i++) | |||
288 | { | |||
289 | /* calc COM */ | |||
290 | clear_rvec(com); | |||
291 | mtot = 0; | |||
292 | for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++) | |||
293 | { | |||
294 | m = atom[j].m; | |||
295 | for (d = 0; d < DIM3; d++) | |||
296 | { | |||
297 | com[d] += m*x[j][d]; | |||
298 | } | |||
299 | mtot += m; | |||
300 | } | |||
301 | /* calculate final COM */ | |||
302 | svmul(1.0/mtot, com, com); | |||
303 | ||||
304 | /* check if COM is outside box */ | |||
305 | copy_rvec(com, new_com); | |||
306 | switch (unitcell_enum) | |||
307 | { | |||
308 | case euRect: | |||
309 | put_atoms_in_box(ePBC, box, 1, &new_com); | |||
310 | break; | |||
311 | case euTric: | |||
312 | put_atoms_in_triclinic_unitcell(ecenter, box, 1, &new_com); | |||
313 | break; | |||
314 | case euCompact: | |||
315 | put_atoms_in_compact_unitcell(ePBC, ecenter, box, 1, &new_com); | |||
316 | break; | |||
317 | } | |||
318 | rvec_sub(new_com, com, shift); | |||
319 | if (norm2(shift) > 0) | |||
320 | { | |||
321 | if (debug) | |||
322 | { | |||
323 | fprintf(debug, "\nShifting position of molecule %d " | |||
324 | "by %8.3f %8.3f %8.3f\n", i+1, | |||
325 | shift[XX0], shift[YY1], shift[ZZ2]); | |||
326 | } | |||
327 | for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++) | |||
328 | { | |||
329 | rvec_inc(x[j], shift); | |||
330 | } | |||
331 | } | |||
332 | } | |||
333 | } | |||
334 | ||||
335 | static void put_residue_com_in_box(int unitcell_enum, int ecenter, | |||
336 | int natoms, t_atom atom[], | |||
337 | int ePBC, matrix box, rvec x[]) | |||
338 | { | |||
339 | atom_id i, j, res_start, res_end, res_nat; | |||
340 | int d, presnr; | |||
341 | real m; | |||
342 | double mtot; | |||
343 | rvec box_center, com, new_com, shift; | |||
344 | ||||
345 | calc_box_center(ecenter, box, box_center); | |||
346 | ||||
347 | presnr = NOTSET-12345; | |||
348 | res_start = 0; | |||
349 | clear_rvec(com); | |||
350 | mtot = 0; | |||
351 | for (i = 0; i < natoms+1; i++) | |||
352 | { | |||
353 | if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET-12345)) | |||
354 | { | |||
355 | /* calculate final COM */ | |||
356 | res_end = i; | |||
357 | res_nat = res_end - res_start; | |||
358 | svmul(1.0/mtot, com, com); | |||
359 | ||||
360 | /* check if COM is outside box */ | |||
361 | copy_rvec(com, new_com); | |||
362 | switch (unitcell_enum) | |||
363 | { | |||
364 | case euRect: | |||
365 | put_atoms_in_box(ePBC, box, 1, &new_com); | |||
366 | break; | |||
367 | case euTric: | |||
368 | put_atoms_in_triclinic_unitcell(ecenter, box, 1, &new_com); | |||
369 | break; | |||
370 | case euCompact: | |||
371 | put_atoms_in_compact_unitcell(ePBC, ecenter, box, 1, &new_com); | |||
372 | break; | |||
373 | } | |||
374 | rvec_sub(new_com, com, shift); | |||
375 | if (norm2(shift)) | |||
376 | { | |||
377 | if (debug) | |||
378 | { | |||
379 | fprintf(debug, "\nShifting position of residue %d (atoms %u-%u) " | |||
380 | "by %g,%g,%g\n", atom[res_start].resind+1, | |||
381 | res_start+1, res_end+1, shift[XX0], shift[YY1], shift[ZZ2]); | |||
382 | } | |||
383 | for (j = res_start; j < res_end; j++) | |||
384 | { | |||
385 | rvec_inc(x[j], shift); | |||
386 | } | |||
387 | } | |||
388 | clear_rvec(com); | |||
389 | mtot = 0; | |||
390 | ||||
391 | /* remember start of new residue */ | |||
392 | res_start = i; | |||
393 | } | |||
394 | if (i < natoms) | |||
395 | { | |||
396 | /* calc COM */ | |||
397 | m = atom[i].m; | |||
398 | for (d = 0; d < DIM3; d++) | |||
399 | { | |||
400 | com[d] += m*x[i][d]; | |||
401 | } | |||
402 | mtot += m; | |||
403 | ||||
404 | presnr = atom[i].resind; | |||
405 | } | |||
406 | } | |||
407 | } | |||
408 | ||||
409 | static void center_x(int ecenter, rvec x[], matrix box, int n, int nc, atom_id ci[]) | |||
410 | { | |||
411 | int i, m, ai; | |||
412 | rvec cmin, cmax, box_center, dx; | |||
413 | ||||
414 | if (nc > 0) | |||
415 | { | |||
416 | copy_rvec(x[ci[0]], cmin); | |||
417 | copy_rvec(x[ci[0]], cmax); | |||
418 | for (i = 0; i < nc; i++) | |||
419 | { | |||
420 | ai = ci[i]; | |||
421 | for (m = 0; m < DIM3; m++) | |||
422 | { | |||
423 | if (x[ai][m] < cmin[m]) | |||
424 | { | |||
425 | cmin[m] = x[ai][m]; | |||
426 | } | |||
427 | else if (x[ai][m] > cmax[m]) | |||
428 | { | |||
429 | cmax[m] = x[ai][m]; | |||
430 | } | |||
431 | } | |||
432 | } | |||
433 | calc_box_center(ecenter, box, box_center); | |||
434 | for (m = 0; m < DIM3; m++) | |||
435 | { | |||
436 | dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5; | |||
437 | } | |||
438 | ||||
439 | for (i = 0; i < n; i++) | |||
440 | { | |||
441 | rvec_inc(x[i], dx); | |||
442 | } | |||
443 | } | |||
444 | } | |||
445 | ||||
446 | static void mk_filenm(char *base, const char *ext, int ndigit, int file_nr, | |||
447 | char out_file[]) | |||
448 | { | |||
449 | char nbuf[128]; | |||
450 | int nd = 0, fnr; | |||
451 | ||||
452 | strcpy(out_file, base); | |||
453 | fnr = file_nr; | |||
454 | do | |||
455 | { | |||
456 | fnr /= 10; | |||
457 | nd++; | |||
458 | } | |||
459 | while (fnr > 0); | |||
460 | ||||
461 | if (nd < ndigit) | |||
462 | { | |||
463 | strncat(out_file, "00000000000", ndigit-nd)__builtin_strncat (out_file, "00000000000", ndigit-nd); | |||
464 | } | |||
465 | sprintf(nbuf, "%d.", file_nr); | |||
466 | strcat(out_file, nbuf); | |||
467 | strcat(out_file, ext); | |||
468 | } | |||
469 | ||||
470 | void check_trn(const char *fn) | |||
471 | { | |||
472 | if ((fn2ftp(fn) != efTRJ) && (fn2ftp(fn) != efTRR)) | |||
473 | { | |||
474 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 474, "%s is not a trajectory file, exiting\n", fn); | |||
475 | } | |||
476 | } | |||
477 | ||||
478 | #ifndef GMX_NATIVE_WINDOWS | |||
479 | void do_trunc(const char *fn, real t0) | |||
480 | { | |||
481 | t_fileio *in; | |||
482 | FILE *fp; | |||
483 | gmx_bool bStop, bOK; | |||
484 | t_trnheader sh; | |||
485 | gmx_off_t fpos; | |||
486 | char yesno[256]; | |||
487 | int j; | |||
488 | real t = 0; | |||
489 | ||||
490 | if (t0 == -1) | |||
491 | { | |||
492 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 492, "You forgot to set the truncation time"); | |||
493 | } | |||
494 | ||||
495 | /* Check whether this is a .trj file */ | |||
496 | check_trn(fn); | |||
497 | ||||
498 | in = open_trn(fn, "r"); | |||
499 | fp = gmx_fio_getfp(in); | |||
500 | if (fp == NULL((void*)0)) | |||
501 | { | |||
502 | fprintf(stderrstderr, "Sorry, can not trunc %s, truncation of this filetype is not supported\n", fn); | |||
503 | close_trn(in); | |||
504 | } | |||
505 | else | |||
506 | { | |||
507 | j = 0; | |||
508 | fpos = gmx_fio_ftell(in); | |||
509 | bStop = FALSE0; | |||
510 | while (!bStop && fread_trnheader(in, &sh, &bOK)) | |||
511 | { | |||
512 | fread_htrn(in, &sh, NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0)); | |||
513 | fpos = gmx_ftell(fp); | |||
514 | t = sh.t; | |||
515 | if (t >= t0) | |||
516 | { | |||
517 | gmx_fseek(fp, fpos, SEEK_SET0); | |||
518 | bStop = TRUE1; | |||
519 | } | |||
520 | } | |||
521 | if (bStop) | |||
522 | { | |||
523 | fprintf(stderrstderr, "Do you REALLY want to truncate this trajectory (%s) at:\n" | |||
524 | "frame %d, time %g, bytes %ld ??? (type YES if so)\n", | |||
525 | fn, j, t, (long int)fpos); | |||
526 | if (1 != scanf("%s", yesno)) | |||
527 | { | |||
528 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 528, "Error reading user input"); | |||
529 | } | |||
530 | if (strcmp(yesno, "YES")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (yesno) && __builtin_constant_p ("YES") && ( __s1_len = strlen (yesno), __s2_len = strlen ("YES"), (!((size_t )(const void *)((yesno) + 1) - (size_t)(const void *)(yesno) == 1) || __s1_len >= 4) && (!((size_t)(const void *) (("YES") + 1) - (size_t)(const void *)("YES") == 1) || __s2_len >= 4)) ? __builtin_strcmp (yesno, "YES") : (__builtin_constant_p (yesno) && ((size_t)(const void *)((yesno) + 1) - (size_t )(const void *)(yesno) == 1) && (__s1_len = strlen (yesno ), __s1_len < 4) ? (__builtin_constant_p ("YES") && ((size_t)(const void *)(("YES") + 1) - (size_t)(const void * )("YES") == 1) ? __builtin_strcmp (yesno, "YES") : (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) ("YES"); int __result = (((const unsigned char *) (const char *) (yesno))[0] - __s2[0]); if (__s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (yesno))[1] - __s2[1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (yesno))[2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (yesno))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p ("YES") && ((size_t)(const void *)(("YES") + 1) - (size_t )(const void *)("YES") == 1) && (__s2_len = strlen ("YES" ), __s2_len < 4) ? (__builtin_constant_p (yesno) && ((size_t)(const void *)((yesno) + 1) - (size_t)(const void * )(yesno) == 1) ? __builtin_strcmp (yesno, "YES") : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (yesno); int __result = (((const unsigned char *) (const char *) ("YES"))[0] - __s2[0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ("YES"))[1] - __s2[1]); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ("YES"))[2] - __s2[2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) ("YES"))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (yesno, "YES")))); }) == 0) | |||
531 | { | |||
532 | fprintf(stderrstderr, "Once again, I'm gonna DO this...\n"); | |||
533 | close_trn(in); | |||
534 | if (0 != truncate(fn, fpos)) | |||
535 | { | |||
536 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 536, "Error truncating file %s", fn); | |||
537 | } | |||
538 | } | |||
539 | else | |||
540 | { | |||
541 | fprintf(stderrstderr, "Ok, I'll forget about it\n"); | |||
542 | } | |||
543 | } | |||
544 | else | |||
545 | { | |||
546 | fprintf(stderrstderr, "Already at end of file (t=%g)...\n", t); | |||
547 | close_trn(in); | |||
548 | } | |||
549 | } | |||
550 | } | |||
551 | #endif | |||
552 | ||||
553 | /*! \brief Read a full molecular topology if useful and available. | |||
554 | * | |||
555 | * If the input trajectory file is not in TNG format, and the output | |||
556 | * file is in TNG format, then we want to try to read a full topology | |||
557 | * (if available), so that we can write molecule information to the | |||
558 | * output file. The full topology provides better molecule information | |||
559 | * than is available from the normal t_topology data used by GROMACS | |||
560 | * tools. | |||
561 | * | |||
562 | * Also, the t_topology is only read under (different) particular | |||
563 | * conditions. If both apply, then a .tpr file might be read | |||
564 | * twice. Trying to fix this redundancy while trjconv is still an | |||
565 | * all-purpose tool does not seem worthwhile. | |||
566 | * | |||
567 | * Because of the way gmx_prepare_tng_writing is implemented, the case | |||
568 | * where the input TNG file has no molecule information will never | |||
569 | * lead to an output TNG file having molecule information. Since | |||
570 | * molecule information will generally be present if the input TNG | |||
571 | * file was written by a GROMACS tool, this seems like reasonable | |||
572 | * behaviour. */ | |||
573 | static gmx_mtop_t *read_mtop_for_tng(const char *tps_file, | |||
574 | const char *input_file, | |||
575 | const char *output_file) | |||
576 | { | |||
577 | gmx_mtop_t *mtop = NULL((void*)0); | |||
578 | ||||
579 | if (fn2bTPX(tps_file) && | |||
580 | efTNG != fn2ftp(input_file) && | |||
581 | efTNG == fn2ftp(output_file)) | |||
582 | { | |||
583 | int temp_natoms = -1; | |||
584 | snew(mtop, 1)(mtop) = save_calloc("mtop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 584, (1), sizeof(*(mtop))); | |||
585 | read_tpx(tps_file, NULL((void*)0), NULL((void*)0), &temp_natoms, | |||
586 | NULL((void*)0), NULL((void*)0), NULL((void*)0), mtop); | |||
587 | } | |||
588 | ||||
589 | return mtop; | |||
590 | } | |||
591 | ||||
592 | int gmx_trjconv(int argc, char *argv[]) | |||
593 | { | |||
594 | const char *desc[] = { | |||
595 | "[THISMODULE] can convert trajectory files in many ways:[BR]", | |||
596 | "* from one format to another[BR]", | |||
597 | "* select a subset of atoms[BR]", | |||
598 | "* change the periodicity representation[BR]", | |||
599 | "* keep multimeric molecules together[BR]", | |||
600 | "* center atoms in the box[BR]", | |||
601 | "* fit atoms to reference structure[BR]", | |||
602 | "* reduce the number of frames[BR]", | |||
603 | "* change the timestamps of the frames ", | |||
604 | "([TT]-t0[tt] and [TT]-timestep[tt])[BR]", | |||
605 | "* cut the trajectory in small subtrajectories according", | |||
606 | "to information in an index file. This allows subsequent analysis of", | |||
607 | "the subtrajectories that could, for example, be the result of a", | |||
608 | "cluster analysis. Use option [TT]-sub[tt].", | |||
609 | "This assumes that the entries in the index file are frame numbers and", | |||
610 | "dumps each group in the index file to a separate trajectory file.[BR]", | |||
611 | "* select frames within a certain range of a quantity given", | |||
612 | "in an [TT].xvg[tt] file.[PAR]", | |||
613 | ||||
614 | "[gmx-trjcat] is better suited for concatenating multiple trajectory files.", | |||
615 | "[PAR]", | |||
616 | ||||
617 | "The following formats are supported for input and output:", | |||
618 | "[TT].xtc[tt], [TT].trr[tt], [TT].trj[tt], [TT].gro[tt], [TT].g96[tt]", | |||
619 | "and [TT].pdb[tt].", | |||
620 | "The file formats are detected from the file extension.", | |||
621 | "The precision of [TT].xtc[tt] and [TT].gro[tt] output is taken from the", | |||
622 | "input file for [TT].xtc[tt], [TT].gro[tt] and [TT].pdb[tt],", | |||
623 | "and from the [TT]-ndec[tt] option for other input formats. The precision", | |||
624 | "is always taken from [TT]-ndec[tt], when this option is set.", | |||
625 | "All other formats have fixed precision. [TT].trr[tt] and [TT].trj[tt]", | |||
626 | "output can be single or double precision, depending on the precision", | |||
627 | "of the [THISMODULE] binary.", | |||
628 | "Note that velocities are only supported in", | |||
629 | "[TT].trr[tt], [TT].trj[tt], [TT].gro[tt] and [TT].g96[tt] files.[PAR]", | |||
630 | ||||
631 | "Option [TT]-sep[tt] can be used to write every frame to a separate", | |||
632 | "[TT].gro, .g96[tt] or [TT].pdb[tt] file. By default, all frames all written to one file.", | |||
633 | "[TT].pdb[tt] files with all frames concatenated can be viewed with", | |||
634 | "[TT]rasmol -nmrpdb[tt].[PAR]", | |||
635 | ||||
636 | "It is possible to select part of your trajectory and write it out", | |||
637 | "to a new trajectory file in order to save disk space, e.g. for leaving", | |||
638 | "out the water from a trajectory of a protein in water.", | |||
639 | "[BB]ALWAYS[bb] put the original trajectory on tape!", | |||
640 | "We recommend to use the portable [TT].xtc[tt] format for your analysis", | |||
641 | "to save disk space and to have portable files.[PAR]", | |||
642 | ||||
643 | "There are two options for fitting the trajectory to a reference", | |||
644 | "either for essential dynamics analysis, etc.", | |||
645 | "The first option is just plain fitting to a reference structure", | |||
646 | "in the structure file. The second option is a progressive fit", | |||
647 | "in which the first timeframe is fitted to the reference structure ", | |||
648 | "in the structure file to obtain and each subsequent timeframe is ", | |||
649 | "fitted to the previously fitted structure. This way a continuous", | |||
650 | "trajectory is generated, which might not be the case when using the", | |||
651 | "regular fit method, e.g. when your protein undergoes large", | |||
652 | "conformational transitions.[PAR]", | |||
653 | ||||
654 | "Option [TT]-pbc[tt] sets the type of periodic boundary condition", | |||
655 | "treatment:[BR]", | |||
656 | "[TT]* mol[tt] puts the center of mass of molecules in the box,", | |||
657 | "and requires a run input file to be supplied with [TT]-s[tt].[BR]", | |||
658 | "[TT]* res[tt] puts the center of mass of residues in the box.[BR]", | |||
659 | "[TT]* atom[tt] puts all the atoms in the box.[BR]", | |||
660 | "[TT]* nojump[tt] checks if atoms jump across the box and then puts", | |||
661 | "them back. This has the effect that all molecules", | |||
662 | "will remain whole (provided they were whole in the initial", | |||
663 | "conformation). [BB]Note[bb] that this ensures a continuous trajectory but", | |||
664 | "molecules may diffuse out of the box. The starting configuration", | |||
665 | "for this procedure is taken from the structure file, if one is", | |||
666 | "supplied, otherwise it is the first frame.[BR]", | |||
667 | "[TT]* cluster[tt] clusters all the atoms in the selected index", | |||
668 | "such that they are all closest to the center of mass of the cluster,", | |||
669 | "which is iteratively updated. [BB]Note[bb] that this will only give meaningful", | |||
670 | "results if you in fact have a cluster. Luckily that can be checked", | |||
671 | "afterwards using a trajectory viewer. Note also that if your molecules", | |||
672 | "are broken this will not work either.[BR]", | |||
673 | "The separate option [TT]-clustercenter[tt] can be used to specify an", | |||
674 | "approximate center for the cluster. This is useful e.g. if you have", | |||
675 | "two big vesicles, and you want to maintain their relative positions.[BR]", | |||
676 | "[TT]* whole[tt] only makes broken molecules whole.[PAR]", | |||
677 | ||||
678 | "Option [TT]-ur[tt] sets the unit cell representation for options", | |||
679 | "[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].", | |||
680 | "All three options give different results for triclinic boxes and", | |||
681 | "identical results for rectangular boxes.", | |||
682 | "[TT]rect[tt] is the ordinary brick shape.", | |||
683 | "[TT]tric[tt] is the triclinic unit cell.", | |||
684 | "[TT]compact[tt] puts all atoms at the closest distance from the center", | |||
685 | "of the box. This can be useful for visualizing e.g. truncated octahedra", | |||
686 | "or rhombic dodecahedra. The center for options [TT]tric[tt] and [TT]compact[tt]", | |||
687 | "is [TT]tric[tt] (see below), unless the option [TT]-boxcenter[tt]", | |||
688 | "is set differently.[PAR]", | |||
689 | ||||
690 | "Option [TT]-center[tt] centers the system in the box. The user can", | |||
691 | "select the group which is used to determine the geometrical center.", | |||
692 | "Option [TT]-boxcenter[tt] sets the location of the center of the box", | |||
693 | "for options [TT]-pbc[tt] and [TT]-center[tt]. The center options are:", | |||
694 | "[TT]tric[tt]: half of the sum of the box vectors,", | |||
695 | "[TT]rect[tt]: half of the box diagonal,", | |||
696 | "[TT]zero[tt]: zero.", | |||
697 | "Use option [TT]-pbc mol[tt] in addition to [TT]-center[tt] when you", | |||
698 | "want all molecules in the box after the centering.[PAR]", | |||
699 | ||||
700 | "It is not always possible to use combinations of [TT]-pbc[tt],", | |||
701 | "[TT]-fit[tt], [TT]-ur[tt] and [TT]-center[tt] to do exactly what", | |||
702 | "you want in one call to [THISMODULE]. Consider using multiple", | |||
703 | "calls, and check out the GROMACS website for suggestions.[PAR]", | |||
704 | ||||
705 | "With [TT]-dt[tt], it is possible to reduce the number of ", | |||
706 | "frames in the output. This option relies on the accuracy of the times", | |||
707 | "in your input trajectory, so if these are inaccurate use the", | |||
708 | "[TT]-timestep[tt] option to modify the time (this can be done", | |||
709 | "simultaneously). For making smooth movies, the program [gmx-filter]", | |||
710 | "can reduce the number of frames while using low-pass frequency", | |||
711 | "filtering, this reduces aliasing of high frequency motions.[PAR]", | |||
712 | ||||
713 | "Using [TT]-trunc[tt] [THISMODULE] can truncate [TT].trj[tt] in place, i.e.", | |||
714 | "without copying the file. This is useful when a run has crashed", | |||
715 | "during disk I/O (i.e. full disk), or when two contiguous", | |||
716 | "trajectories must be concatenated without having double frames.[PAR]", | |||
717 | ||||
718 | "Option [TT]-dump[tt] can be used to extract a frame at or near", | |||
719 | "one specific time from your trajectory.[PAR]", | |||
720 | ||||
721 | "Option [TT]-drop[tt] reads an [TT].xvg[tt] file with times and values.", | |||
722 | "When options [TT]-dropunder[tt] and/or [TT]-dropover[tt] are set,", | |||
723 | "frames with a value below and above the value of the respective options", | |||
724 | "will not be written." | |||
725 | }; | |||
726 | ||||
727 | int pbc_enum; | |||
728 | enum | |||
729 | { | |||
730 | epSel, | |||
731 | epNone, | |||
732 | epComMol, | |||
733 | epComRes, | |||
734 | epComAtom, | |||
735 | epNojump, | |||
736 | epCluster, | |||
737 | epWhole, | |||
738 | epNR | |||
739 | }; | |||
740 | const char *pbc_opt[epNR + 1] = | |||
741 | { | |||
742 | NULL((void*)0), "none", "mol", "res", "atom", "nojump", "cluster", "whole", | |||
743 | NULL((void*)0) | |||
744 | }; | |||
745 | ||||
746 | int unitcell_enum; | |||
747 | const char *unitcell_opt[euNR+1] = | |||
748 | { NULL((void*)0), "rect", "tric", "compact", NULL((void*)0) }; | |||
749 | ||||
750 | enum | |||
751 | { | |||
752 | ecSel, ecTric, ecRect, ecZero, ecNR | |||
753 | }; | |||
754 | const char *center_opt[ecNR+1] = | |||
755 | { NULL((void*)0), "tric", "rect", "zero", NULL((void*)0) }; | |||
756 | int ecenter; | |||
757 | ||||
758 | int fit_enum; | |||
759 | enum | |||
760 | { | |||
761 | efSel, efNone, efFit, efFitXY, efReset, efResetXY, efPFit, efNR | |||
762 | }; | |||
763 | const char *fit[efNR + 1] = | |||
764 | { | |||
765 | NULL((void*)0), "none", "rot+trans", "rotxy+transxy", "translation", "transxy", | |||
766 | "progressive", NULL((void*)0) | |||
767 | }; | |||
768 | ||||
769 | static gmx_bool bSeparate = FALSE0, bVels = TRUE1, bForce = FALSE0, bCONECT = FALSE0; | |||
770 | static gmx_bool bCenter = FALSE0; | |||
771 | static int skip_nr = 1, ndec = 3, nzero = 0; | |||
772 | static real tzero = 0, delta_t = 0, timestep = 0, ttrunc = -1, tdump = -1, split_t = 0; | |||
773 | static rvec newbox = {0, 0, 0}, shift = {0, 0, 0}, trans = {0, 0, 0}; | |||
774 | static char *exec_command = NULL((void*)0); | |||
775 | static real dropunder = 0, dropover = 0; | |||
776 | static gmx_bool bRound = FALSE0; | |||
777 | ||||
778 | t_pargs | |||
779 | pa[] = | |||
780 | { | |||
781 | { "-skip", FALSE0, etINT, | |||
782 | { &skip_nr }, "Only write every nr-th frame" }, | |||
783 | { "-dt", FALSE0, etTIME, | |||
784 | { &delta_t }, | |||
785 | "Only write frame when t MOD dt = first time (%t)" }, | |||
786 | { "-round", FALSE0, etBOOL, | |||
787 | { &bRound }, "Round measurements to nearest picosecond"}, | |||
788 | { "-dump", FALSE0, etTIME, | |||
789 | { &tdump }, "Dump frame nearest specified time (%t)" }, | |||
790 | { "-t0", FALSE0, etTIME, | |||
791 | { &tzero }, | |||
792 | "Starting time (%t) (default: don't change)" }, | |||
793 | { "-timestep", FALSE0, etTIME, | |||
794 | { ×tep }, | |||
795 | "Change time step between input frames (%t)" }, | |||
796 | { "-pbc", FALSE0, etENUM, | |||
797 | { pbc_opt }, | |||
798 | "PBC treatment (see help text for full description)" }, | |||
799 | { "-ur", FALSE0, etENUM, | |||
800 | { unitcell_opt }, "Unit-cell representation" }, | |||
801 | { "-center", FALSE0, etBOOL, | |||
802 | { &bCenter }, "Center atoms in box" }, | |||
803 | { "-boxcenter", FALSE0, etENUM, | |||
804 | { center_opt }, "Center for -pbc and -center" }, | |||
805 | { "-box", FALSE0, etRVEC, | |||
806 | { newbox }, | |||
807 | "Size for new cubic box (default: read from input)" }, | |||
808 | { "-trans", FALSE0, etRVEC, | |||
809 | { trans }, | |||
810 | "All coordinates will be translated by trans. This " | |||
811 | "can advantageously be combined with -pbc mol -ur " | |||
812 | "compact." }, | |||
813 | { "-shift", FALSE0, etRVEC, | |||
814 | { shift }, | |||
815 | "All coordinates will be shifted by framenr*shift" }, | |||
816 | { "-fit", FALSE0, etENUM, | |||
817 | { fit }, | |||
818 | "Fit molecule to ref structure in the structure file" }, | |||
819 | { "-ndec", FALSE0, etINT, | |||
820 | { &ndec }, | |||
821 | "Precision for .xtc and .gro writing in number of " | |||
822 | "decimal places" }, | |||
823 | { "-vel", FALSE0, etBOOL, | |||
824 | { &bVels }, "Read and write velocities if possible" }, | |||
825 | { "-force", FALSE0, etBOOL, | |||
826 | { &bForce }, "Read and write forces if possible" }, | |||
827 | #ifndef GMX_NATIVE_WINDOWS | |||
828 | { "-trunc", FALSE0, etTIME, | |||
829 | { &ttrunc }, | |||
830 | "Truncate input trajectory file after this time (%t)" }, | |||
831 | #endif | |||
832 | { "-exec", FALSE0, etSTR, | |||
833 | { &exec_command }, | |||
834 | "Execute command for every output frame with the " | |||
835 | "frame number as argument" }, | |||
836 | { "-split", FALSE0, etTIME, | |||
837 | { &split_t }, | |||
838 | "Start writing new file when t MOD split = first " | |||
839 | "time (%t)" }, | |||
840 | { "-sep", FALSE0, etBOOL, | |||
841 | { &bSeparate }, | |||
842 | "Write each frame to a separate .gro, .g96 or .pdb " | |||
843 | "file" }, | |||
844 | { "-nzero", FALSE0, etINT, | |||
845 | { &nzero }, | |||
846 | "If the -sep flag is set, use these many digits " | |||
847 | "for the file numbers and prepend zeros as needed" }, | |||
848 | { "-dropunder", FALSE0, etREAL, | |||
849 | { &dropunder }, "Drop all frames below this value" }, | |||
850 | { "-dropover", FALSE0, etREAL, | |||
851 | { &dropover }, "Drop all frames above this value" }, | |||
852 | { "-conect", FALSE0, etBOOL, | |||
853 | { &bCONECT }, | |||
854 | "Add conect records when writing [TT].pdb[tt] files. Useful " | |||
855 | "for visualization of non-standard molecules, e.g. " | |||
856 | "coarse grained ones" } | |||
857 | }; | |||
858 | #define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))) | |||
859 | ||||
860 | FILE *out = NULL((void*)0); | |||
861 | t_trxstatus *trxout = NULL((void*)0); | |||
862 | t_trxstatus *trxin; | |||
863 | int ftp, ftpin = 0, file_nr; | |||
864 | t_trxframe fr, frout; | |||
865 | int flags; | |||
866 | rvec *xmem = NULL((void*)0), *vmem = NULL((void*)0), *fmem = NULL((void*)0); | |||
867 | rvec *xp = NULL((void*)0), x_shift, hbox, box_center, dx; | |||
868 | real xtcpr, lambda, *w_rls = NULL((void*)0); | |||
869 | int m, i, d, frame, outframe, natoms, nout, ncent, nre, newstep = 0, model_nr; | |||
870 | #define SKIP10 10 | |||
871 | t_topology top; | |||
872 | gmx_mtop_t *mtop = NULL((void*)0); | |||
873 | gmx_conect gc = NULL((void*)0); | |||
874 | int ePBC = -1; | |||
875 | t_atoms *atoms = NULL((void*)0), useatoms; | |||
876 | matrix top_box; | |||
877 | atom_id *index, *cindex; | |||
878 | char *grpnm; | |||
| ||||
879 | int *frindex, nrfri; | |||
880 | char *frname; | |||
881 | int ifit, irms, my_clust = -1; | |||
882 | atom_id *ind_fit, *ind_rms; | |||
883 | char *gn_fit, *gn_rms; | |||
884 | t_cluster_ndx *clust = NULL((void*)0); | |||
885 | t_trxstatus **clust_status = NULL((void*)0); | |||
886 | int *clust_status_id = NULL((void*)0); | |||
887 | int ntrxopen = 0; | |||
888 | int *nfwritten = NULL((void*)0); | |||
889 | int ndrop = 0, ncol, drop0 = 0, drop1 = 0, dropuse = 0; | |||
890 | double **dropval; | |||
891 | real tshift = 0, t0 = -1, dt = 0.001, prec; | |||
892 | gmx_bool bFit, bFitXY, bPFit, bReset; | |||
893 | int nfitdim; | |||
894 | gmx_rmpbc_t gpbc = NULL((void*)0); | |||
895 | gmx_bool bRmPBC, bPBCWhole, bPBCcomRes, bPBCcomMol, bPBCcomAtom, bPBC, bNoJump, bCluster; | |||
896 | gmx_bool bCopy, bDoIt, bIndex, bTDump, bSetTime, bTPS = FALSE0, bDTset = FALSE0; | |||
897 | gmx_bool bExec, bTimeStep = FALSE0, bDumpFrame = FALSE0, bSetPrec, bNeedPrec; | |||
898 | gmx_bool bHaveFirstFrame, bHaveNextFrame, bSetBox, bSetUR, bSplit = FALSE0; | |||
899 | gmx_bool bSubTraj = FALSE0, bDropUnder = FALSE0, bDropOver = FALSE0, bTrans = FALSE0; | |||
900 | gmx_bool bWriteFrame, bSplitHere; | |||
901 | const char *top_file, *in_file, *out_file = NULL((void*)0); | |||
902 | char out_file2[256], *charpt; | |||
903 | char *outf_base = NULL((void*)0); | |||
904 | const char *outf_ext = NULL((void*)0); | |||
905 | char top_title[256], title[256], command[256], filemode[5]; | |||
906 | int xdr = 0; | |||
907 | gmx_bool bWarnCompact = FALSE0; | |||
908 | const char *warn; | |||
909 | output_env_t oenv; | |||
910 | ||||
911 | t_filenm fnm[] = { | |||
912 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, | |||
913 | { efTRO, "-o", NULL((void*)0), ffWRITE1<<2 }, | |||
914 | { efTPS, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
915 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
916 | { efNDX, "-fr", "frames", ffOPTRD(1<<1 | 1<<3) }, | |||
917 | { efNDX, "-sub", "cluster", ffOPTRD(1<<1 | 1<<3) }, | |||
918 | { efXVG, "-drop", "drop", ffOPTRD(1<<1 | 1<<3) } | |||
919 | }; | |||
920 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) | |||
921 | ||||
922 | if (!parse_common_args(&argc, argv, | |||
923 | PCA_CAN_BEGIN(1<<6) | PCA_CAN_END(1<<7) | PCA_CAN_VIEW(1<<5) | | |||
924 | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13), | |||
925 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, | |||
926 | 0, NULL((void*)0), &oenv)) | |||
927 | { | |||
928 | return 0; | |||
929 | } | |||
930 | ||||
931 | top_file = ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
932 | init_top(&top); | |||
933 | ||||
934 | /* Check command line */ | |||
935 | in_file = opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
936 | ||||
937 | if (ttrunc != -1) | |||
938 | { | |||
939 | #ifndef GMX_NATIVE_WINDOWS | |||
940 | do_trunc(in_file, ttrunc); | |||
941 | #endif | |||
942 | } | |||
943 | else | |||
944 | { | |||
945 | /* mark active cmdline options */ | |||
946 | bSetBox = opt2parg_bSet("-box", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
947 | bSetTime = opt2parg_bSet("-t0", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
948 | bSetPrec = opt2parg_bSet("-ndec", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
949 | bSetUR = opt2parg_bSet("-ur", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
950 | bExec = opt2parg_bSet("-exec", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
951 | bTimeStep = opt2parg_bSet("-timestep", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
952 | bTDump = opt2parg_bSet("-dump", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
953 | bDropUnder = opt2parg_bSet("-dropunder", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
954 | bDropOver = opt2parg_bSet("-dropover", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
955 | bTrans = opt2parg_bSet("-trans", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); | |||
956 | bSplit = (split_t != 0); | |||
957 | ||||
958 | /* parse enum options */ | |||
959 | fit_enum = nenum(fit); | |||
960 | bFit = (fit_enum == efFit || fit_enum == efFitXY); | |||
961 | bFitXY = fit_enum == efFitXY; | |||
962 | bReset = (fit_enum == efReset || fit_enum == efResetXY); | |||
963 | bPFit = fit_enum == efPFit; | |||
964 | pbc_enum = nenum(pbc_opt); | |||
965 | bPBCWhole = pbc_enum == epWhole; | |||
966 | bPBCcomRes = pbc_enum == epComRes; | |||
967 | bPBCcomMol = pbc_enum == epComMol; | |||
968 | bPBCcomAtom = pbc_enum == epComAtom; | |||
969 | bNoJump = pbc_enum == epNojump; | |||
970 | bCluster = pbc_enum == epCluster; | |||
971 | bPBC = pbc_enum != epNone; | |||
972 | unitcell_enum = nenum(unitcell_opt); | |||
973 | ecenter = nenum(center_opt) - ecTric; | |||
974 | ||||
975 | /* set and check option dependencies */ | |||
976 | if (bPFit) | |||
977 | { | |||
978 | bFit = TRUE1; /* for pfit, fit *must* be set */ | |||
979 | } | |||
980 | if (bFit) | |||
981 | { | |||
982 | bReset = TRUE1; /* for fit, reset *must* be set */ | |||
983 | } | |||
984 | nfitdim = 0; | |||
985 | if (bFit || bReset) | |||
986 | { | |||
987 | nfitdim = (fit_enum == efFitXY || fit_enum == efResetXY) ? 2 : 3; | |||
988 | } | |||
989 | bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol; | |||
990 | ||||
991 | if (bSetUR) | |||
992 | { | |||
993 | if (!(bPBCcomRes || bPBCcomMol || bPBCcomAtom)) | |||
994 | { | |||
995 | fprintf(stderrstderr, | |||
996 | "WARNING: Option for unitcell representation (-ur %s)\n" | |||
997 | " only has effect in combination with -pbc %s, %s or %s.\n" | |||
998 | " Ingoring unitcell representation.\n\n", | |||
999 | unitcell_opt[0], pbc_opt[2], pbc_opt[3], pbc_opt[4]); | |||
1000 | bSetUR = FALSE0; | |||
1001 | } | |||
1002 | } | |||
1003 | if (bFit && bPBC) | |||
1004 | { | |||
1005 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1005, "PBC condition treatment does not work together with rotational fit.\n" | |||
1006 | "Please do the PBC condition treatment first and then run trjconv in a second step\n" | |||
1007 | "for the rotational fit.\n" | |||
1008 | "First doing the rotational fit and then doing the PBC treatment gives incorrect\n" | |||
1009 | "results!"); | |||
1010 | } | |||
1011 | ||||
1012 | /* ndec is in nr of decimal places, prec is a multiplication factor: */ | |||
1013 | prec = 1; | |||
1014 | for (i = 0; i < ndec; i++) | |||
1015 | { | |||
1016 | prec *= 10; | |||
1017 | } | |||
1018 | ||||
1019 | bIndex = ftp2bSet(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
1020 | ||||
1021 | ||||
1022 | /* Determine output type */ | |||
1023 | out_file = opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
1024 | ftp = fn2ftp(out_file); | |||
1025 | fprintf(stderrstderr, "Will write %s: %s\n", ftp2ext(ftp), ftp2desc(ftp)); | |||
1026 | bNeedPrec = (ftp == efXTC || ftp == efGRO); | |||
1027 | if (bVels) | |||
1028 | { | |||
1029 | /* check if velocities are possible in input and output files */ | |||
1030 | ftpin = fn2ftp(in_file); | |||
1031 | bVels = (ftp == efTRR || ftp == efTRJ || ftp == efGRO || ftp == efG96) | |||
1032 | && (ftpin == efTRR || ftpin == efTRJ || ftpin == efGRO || ftpin == efG96 || | |||
1033 | ftpin == efCPT); | |||
1034 | } | |||
1035 | if (bSeparate || bSplit) | |||
1036 | { | |||
1037 | outf_ext = strrchr(out_file, '.'); | |||
1038 | if (outf_ext == NULL((void*)0)) | |||
1039 | { | |||
1040 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1040, "Output file name '%s' does not contain a '.'", out_file); | |||
1041 | } | |||
1042 | outf_base = strdup(out_file)(__extension__ (__builtin_constant_p (out_file) && (( size_t)(const void *)((out_file) + 1) - (size_t)(const void * )(out_file) == 1) ? (((const char *) (out_file))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (out_file) + 1; char *__retval = (char *) malloc (__len ); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval , out_file, __len); __retval; })) : __strdup (out_file))); | |||
1043 | outf_base[outf_ext - out_file] = '\0'; | |||
1044 | } | |||
1045 | ||||
1046 | bSubTraj = opt2bSet("-sub", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
1047 | if (bSubTraj) | |||
1048 | { | |||
1049 | if ((ftp != efXTC) && (ftp != efTRR)) | |||
1050 | { | |||
1051 | /* It seems likely that other trajectory file types | |||
1052 | * could work here. */ | |||
1053 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1053, "Can only use the sub option with output file types " | |||
1054 | "xtc and trr"); | |||
1055 | } | |||
1056 | clust = cluster_index(NULL((void*)0), opt2fn("-sub", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); | |||
1057 | ||||
1058 | /* Check for number of files disabled, as FOPEN_MAX is not the correct | |||
1059 | * number to check for. In my linux box it is only 16. | |||
1060 | */ | |||
1061 | if (0 && (clust->clust->nr > FOPEN_MAX16-4)) | |||
1062 | { | |||
1063 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1063, "Can not open enough (%d) files to write all the" | |||
1064 | " trajectories.\ntry splitting the index file in %d parts.\n" | |||
1065 | "FOPEN_MAX = %d", | |||
1066 | clust->clust->nr, 1+clust->clust->nr/FOPEN_MAX16, FOPEN_MAX16); | |||
1067 | } | |||
1068 | gmx_warning("The -sub option could require as many open output files as there are\n" | |||
1069 | "index groups in the file (%d). If you get I/O errors opening new files,\n" | |||
1070 | "try reducing the number of index groups in the file, and perhaps\n" | |||
1071 | "using trjconv -sub several times on different chunks of your index file.\n", | |||
1072 | clust->clust->nr); | |||
1073 | ||||
1074 | snew(clust_status, clust->clust->nr)(clust_status) = save_calloc("clust_status", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1074, (clust->clust->nr), sizeof(*(clust_status))); | |||
1075 | snew(clust_status_id, clust->clust->nr)(clust_status_id) = save_calloc("clust_status_id", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1075, (clust->clust->nr), sizeof(*(clust_status_id))); | |||
1076 | snew(nfwritten, clust->clust->nr)(nfwritten) = save_calloc("nfwritten", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1076, (clust->clust->nr), sizeof(*(nfwritten))); | |||
1077 | for (i = 0; (i < clust->clust->nr); i++) | |||
1078 | { | |||
1079 | clust_status[i] = NULL((void*)0); | |||
1080 | clust_status_id[i] = -1; | |||
1081 | } | |||
1082 | bSeparate = bSplit = FALSE0; | |||
1083 | } | |||
1084 | /* skipping */ | |||
1085 | if (skip_nr <= 0) | |||
1086 | { | |||
1087 | } | |||
1088 | ||||
1089 | mtop = read_mtop_for_tng(top_file, in_file, out_file); | |||
1090 | ||||
1091 | /* Determine whether to read a topology */ | |||
1092 | bTPS = (ftp2bSet(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || | |||
1093 | bRmPBC || bReset || bPBCcomMol || bCluster || | |||
1094 | (ftp == efGRO) || (ftp == efPDB) || bCONECT); | |||
1095 | ||||
1096 | /* Determine if when can read index groups */ | |||
1097 | bIndex = (bIndex || bTPS); | |||
1098 | ||||
1099 | if (bTPS) | |||
1100 | { | |||
1101 | read_tps_conf(top_file, top_title, &top, &ePBC, &xp, NULL((void*)0), top_box, | |||
1102 | bReset || bPBCcomRes); | |||
1103 | atoms = &top.atoms; | |||
1104 | ||||
1105 | if (0 == top.mols.nr && (bCluster || bPBCcomMol)) | |||
1106 | { | |||
1107 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1107, "Option -pbc %s requires a .tpr file for the -s option", pbc_opt[pbc_enum]); | |||
1108 | } | |||
1109 | ||||
1110 | /* top_title is only used for gro and pdb, | |||
1111 | * the header in such a file is top_title t= ... | |||
1112 | * to prevent a double t=, remove it from top_title | |||
1113 | */ | |||
1114 | if ((charpt = strstr(top_title, " t= "))) | |||
1115 | { | |||
1116 | charpt[0] = '\0'; | |||
1117 | } | |||
1118 | ||||
1119 | if (bCONECT) | |||
1120 | { | |||
1121 | gc = gmx_conect_generate(&top); | |||
1122 | } | |||
1123 | if (bRmPBC) | |||
1124 | { | |||
1125 | gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr); | |||
1126 | } | |||
1127 | } | |||
1128 | ||||
1129 | /* get frame number index */ | |||
1130 | frindex = NULL((void*)0); | |||
1131 | if (opt2bSet("-fr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) | |||
1132 | { | |||
1133 | printf("Select groups of frame number indices:\n"); | |||
1134 | rd_index(opt2fn("-fr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &nrfri, (atom_id **)&frindex, &frname); | |||
1135 | if (debug) | |||
1136 | { | |||
1137 | for (i = 0; i < nrfri; i++) | |||
1138 | { | |||
1139 | fprintf(debug, "frindex[%4d]=%4d\n", i, frindex[i]); | |||
1140 | } | |||
1141 | } | |||
1142 | } | |||
1143 | ||||
1144 | /* get index groups etc. */ | |||
1145 | if (bReset) | |||
1146 | { | |||
1147 | printf("Select group for %s fit\n", | |||
1148 | bFit ? "least squares" : "translational"); | |||
1149 | get_index(atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
1150 | 1, &ifit, &ind_fit, &gn_fit); | |||
1151 | ||||
1152 | if (bFit) | |||
1153 | { | |||
1154 | if (ifit < 2) | |||
1155 | { | |||
1156 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1156, "Need at least 2 atoms to fit!\n"); | |||
1157 | } | |||
1158 | else if (ifit == 3) | |||
1159 | { | |||
1160 | fprintf(stderrstderr, "WARNING: fitting with only 2 atoms is not unique\n"); | |||
1161 | } | |||
1162 | } | |||
1163 | } | |||
1164 | else if (bCluster) | |||
1165 | { | |||
1166 | printf("Select group for clustering\n"); | |||
1167 | get_index(atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
1168 | 1, &ifit, &ind_fit, &gn_fit); | |||
1169 | } | |||
1170 | ||||
1171 | if (bIndex) | |||
1172 | { | |||
1173 | if (bCenter) | |||
1174 | { | |||
1175 | printf("Select group for centering\n"); | |||
1176 | get_index(atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
1177 | 1, &ncent, &cindex, &grpnm); | |||
1178 | } | |||
1179 | printf("Select group for output\n"); | |||
1180 | get_index(atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
1181 | 1, &nout, &index, &grpnm); | |||
1182 | } | |||
1183 | else | |||
1184 | { | |||
1185 | /* no index file, so read natoms from TRX */ | |||
1186 | if (!read_first_frame(oenv, &trxin, in_file, &fr, TRX_DONT_SKIP(1<<6))) | |||
1187 | { | |||
1188 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1188, "Could not read a frame from %s", in_file); | |||
1189 | } | |||
1190 | natoms = fr.natoms; | |||
1191 | close_trj(trxin); | |||
1192 | sfree(fr.x)save_free("fr.x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1192, (fr.x)); | |||
1193 | snew(index, natoms)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1193, (natoms), sizeof(*(index))); | |||
1194 | for (i = 0; i < natoms; i++) | |||
1195 | { | |||
1196 | index[i] = i; | |||
1197 | } | |||
1198 | nout = natoms; | |||
1199 | if (bCenter) | |||
1200 | { | |||
1201 | ncent = nout; | |||
1202 | cindex = index; | |||
1203 | } | |||
1204 | } | |||
1205 | ||||
1206 | if (bReset) | |||
1207 | { | |||
1208 | snew(w_rls, atoms->nr)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1208, (atoms->nr), sizeof(*(w_rls))); | |||
1209 | for (i = 0; (i < ifit); i++) | |||
1210 | { | |||
1211 | w_rls[ind_fit[i]] = atoms->atom[ind_fit[i]].m; | |||
1212 | } | |||
1213 | ||||
1214 | /* Restore reference structure and set to origin, | |||
1215 | store original location (to put structure back) */ | |||
1216 | if (bRmPBC) | |||
1217 | { | |||
1218 | gmx_rmpbc(gpbc, top.atoms.nr, top_box, xp); | |||
1219 | } | |||
1220 | copy_rvec(xp[index[0]], x_shift); | |||
1221 | reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, NULL((void*)0), xp, w_rls); | |||
1222 | rvec_dec(x_shift, xp[index[0]]); | |||
1223 | } | |||
1224 | else | |||
1225 | { | |||
1226 | clear_rvec(x_shift); | |||
1227 | } | |||
1228 | ||||
1229 | if (bDropUnder || bDropOver) | |||
1230 | { | |||
1231 | /* Read the .xvg file with the drop values */ | |||
1232 | fprintf(stderrstderr, "\nReading drop file ..."); | |||
1233 | ndrop = read_xvg(opt2fn("-drop", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &dropval, &ncol); | |||
1234 | fprintf(stderrstderr, " %d time points\n", ndrop); | |||
1235 | if (ndrop == 0 || ncol < 2) | |||
1236 | { | |||
1237 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1237, "Found no data points in %s", | |||
1238 | opt2fn("-drop", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); | |||
1239 | } | |||
1240 | drop0 = 0; | |||
1241 | drop1 = 0; | |||
1242 | } | |||
1243 | ||||
1244 | /* Make atoms struct for output in GRO or PDB files */ | |||
1245 | if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB)) | |||
1246 | { | |||
1247 | /* get memory for stuff to go in .pdb file */ | |||
1248 | init_t_atoms(&useatoms, atoms->nr, FALSE0); | |||
1249 | sfree(useatoms.resinfo)save_free("useatoms.resinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1249, (useatoms.resinfo)); | |||
1250 | useatoms.resinfo = atoms->resinfo; | |||
1251 | for (i = 0; (i < nout); i++) | |||
1252 | { | |||
1253 | useatoms.atomname[i] = atoms->atomname[index[i]]; | |||
1254 | useatoms.atom[i] = atoms->atom[index[i]]; | |||
1255 | useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1)(((useatoms.nres) > (useatoms.atom[i].resind+1)) ? (useatoms .nres) : (useatoms.atom[i].resind+1) ); | |||
1256 | } | |||
1257 | useatoms.nr = nout; | |||
1258 | } | |||
1259 | /* select what to read */ | |||
1260 | if (ftp == efTRR || ftp == efTRJ) | |||
1261 | { | |||
1262 | flags = TRX_READ_X(1<<0); | |||
1263 | } | |||
1264 | else | |||
1265 | { | |||
1266 | flags = TRX_NEED_X(1<<1); | |||
1267 | } | |||
1268 | if (bVels) | |||
1269 | { | |||
1270 | flags = flags | TRX_READ_V(1<<2); | |||
1271 | } | |||
1272 | if (bForce) | |||
1273 | { | |||
1274 | flags = flags | TRX_READ_F(1<<4); | |||
1275 | } | |||
1276 | ||||
1277 | /* open trx file for reading */ | |||
1278 | bHaveFirstFrame = read_first_frame(oenv, &trxin, in_file, &fr, flags); | |||
1279 | if (fr.bPrec) | |||
1280 | { | |||
1281 | fprintf(stderrstderr, "\nPrecision of %s is %g (nm)\n", in_file, 1/fr.prec); | |||
1282 | } | |||
1283 | if (bNeedPrec) | |||
1284 | { | |||
1285 | if (bSetPrec || !fr.bPrec) | |||
1286 | { | |||
1287 | fprintf(stderrstderr, "\nSetting output precision to %g (nm)\n", 1/prec); | |||
1288 | } | |||
1289 | else | |||
1290 | { | |||
1291 | fprintf(stderrstderr, "Using output precision of %g (nm)\n", 1/prec); | |||
1292 | } | |||
1293 | } | |||
1294 | ||||
1295 | if (bHaveFirstFrame) | |||
1296 | { | |||
1297 | set_trxframe_ePBC(&fr, ePBC); | |||
1298 | ||||
1299 | natoms = fr.natoms; | |||
1300 | ||||
1301 | if (bSetTime) | |||
1302 | { | |||
1303 | tshift = tzero-fr.time; | |||
1304 | } | |||
1305 | else | |||
1306 | { | |||
1307 | tzero = fr.time; | |||
1308 | } | |||
1309 | ||||
1310 | bCopy = FALSE0; | |||
1311 | if (bIndex) | |||
1312 | { | |||
1313 | /* check if index is meaningful */ | |||
1314 | for (i = 0; i < nout; i++) | |||
1315 | { | |||
1316 | if (index[i] >= natoms) | |||
1317 | { | |||
1318 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1318, | |||
1319 | "Index[%d] %d is larger than the number of atoms in the\n" | |||
1320 | "trajectory file (%d). There is a mismatch in the contents\n" | |||
1321 | "of your -f, -s and/or -n files.", i, index[i]+1, natoms); | |||
1322 | } | |||
1323 | bCopy = bCopy || (i != index[i]); | |||
1324 | } | |||
1325 | } | |||
1326 | ||||
1327 | /* open output for writing */ | |||
1328 | strcpy(filemode, "w"); | |||
1329 | switch (ftp) | |||
1330 | { | |||
1331 | case efTNG: | |||
1332 | trjtools_gmx_prepare_tng_writing(out_file, | |||
| ||||
1333 | filemode[0], | |||
1334 | trxin, | |||
1335 | &trxout, | |||
1336 | NULL((void*)0), | |||
1337 | nout, | |||
1338 | mtop, | |||
1339 | index, | |||
1340 | grpnm); | |||
1341 | break; | |||
1342 | case efXTC: | |||
1343 | case efTRR: | |||
1344 | case efTRJ: | |||
1345 | out = NULL((void*)0); | |||
1346 | if (!bSplit && !bSubTraj) | |||
1347 | { | |||
1348 | trxout = open_trx(out_file, filemode); | |||
1349 | } | |||
1350 | break; | |||
1351 | case efGRO: | |||
1352 | case efG96: | |||
1353 | case efPDB: | |||
1354 | if (( !bSeparate && !bSplit ) && !bSubTraj) | |||
1355 | { | |||
1356 | out = gmx_ffopen(out_file, filemode); | |||
1357 | } | |||
1358 | break; | |||
1359 | default: | |||
1360 | gmx_incons("Illegal output file format")_gmx_error("incons", "Illegal output file format", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1360); | |||
1361 | } | |||
1362 | ||||
1363 | if (bCopy) | |||
1364 | { | |||
1365 | snew(xmem, nout)(xmem) = save_calloc("xmem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1365, (nout), sizeof(*(xmem))); | |||
1366 | if (bVels) | |||
1367 | { | |||
1368 | snew(vmem, nout)(vmem) = save_calloc("vmem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1368, (nout), sizeof(*(vmem))); | |||
1369 | } | |||
1370 | if (bForce) | |||
1371 | { | |||
1372 | snew(fmem, nout)(fmem) = save_calloc("fmem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1372, (nout), sizeof(*(fmem))); | |||
1373 | } | |||
1374 | } | |||
1375 | ||||
1376 | /* Start the big loop over frames */ | |||
1377 | file_nr = 0; | |||
1378 | frame = 0; | |||
1379 | outframe = 0; | |||
1380 | model_nr = 0; | |||
1381 | bDTset = FALSE0; | |||
1382 | ||||
1383 | /* Main loop over frames */ | |||
1384 | do | |||
1385 | { | |||
1386 | if (!fr.bStep) | |||
1387 | { | |||
1388 | /* set the step */ | |||
1389 | fr.step = newstep; | |||
1390 | newstep++; | |||
1391 | } | |||
1392 | if (bSubTraj) | |||
1393 | { | |||
1394 | /*if (frame >= clust->clust->nra) | |||
1395 | gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file\n");*/ | |||
1396 | if (frame > clust->maxframe) | |||
1397 | { | |||
1398 | my_clust = -1; | |||
1399 | } | |||
1400 | else | |||
1401 | { | |||
1402 | my_clust = clust->inv_clust[frame]; | |||
1403 | } | |||
1404 | if ((my_clust < 0) || (my_clust >= clust->clust->nr) || | |||
1405 | (my_clust == NO_ATID(atom_id)(~0))) | |||
1406 | { | |||
1407 | my_clust = -1; | |||
1408 | } | |||
1409 | } | |||
1410 | ||||
1411 | if (bSetBox) | |||
1412 | { | |||
1413 | /* generate new box */ | |||
1414 | clear_mat(fr.box); | |||
1415 | for (m = 0; m < DIM3; m++) | |||
1416 | { | |||
1417 | fr.box[m][m] = newbox[m]; | |||
1418 | } | |||
1419 | } | |||
1420 | ||||
1421 | if (bTrans) | |||
1422 | { | |||
1423 | for (i = 0; i < natoms; i++) | |||
1424 | { | |||
1425 | rvec_inc(fr.x[i], trans); | |||
1426 | } | |||
1427 | } | |||
1428 | ||||
1429 | if (bTDump) | |||
1430 | { | |||
1431 | /* determine timestep */ | |||
1432 | if (t0 == -1) | |||
1433 | { | |||
1434 | t0 = fr.time; | |||
1435 | } | |||
1436 | else | |||
1437 | { | |||
1438 | if (!bDTset) | |||
1439 | { | |||
1440 | dt = fr.time-t0; | |||
1441 | bDTset = TRUE1; | |||
1442 | } | |||
1443 | } | |||
1444 | /* This is not very elegant, as one can not dump a frame after | |||
1445 | * a timestep with is more than twice as small as the first one. */ | |||
1446 | bDumpFrame = (fr.time > tdump-0.5*dt) && (fr.time <= tdump+0.5*dt); | |||
1447 | } | |||
1448 | else | |||
1449 | { | |||
1450 | bDumpFrame = FALSE0; | |||
1451 | } | |||
1452 | ||||
1453 | /* determine if an atom jumped across the box and reset it if so */ | |||
1454 | if (bNoJump && (bTPS || frame != 0)) | |||
1455 | { | |||
1456 | for (d = 0; d < DIM3; d++) | |||
1457 | { | |||
1458 | hbox[d] = 0.5*fr.box[d][d]; | |||
1459 | } | |||
1460 | for (i = 0; i < natoms; i++) | |||
1461 | { | |||
1462 | if (bReset) | |||
1463 | { | |||
1464 | rvec_dec(fr.x[i], x_shift); | |||
1465 | } | |||
1466 | for (m = DIM3-1; m >= 0; m--) | |||
1467 | { | |||
1468 | if (hbox[m] > 0) | |||
1469 | { | |||
1470 | while (fr.x[i][m]-xp[i][m] <= -hbox[m]) | |||
1471 | { | |||
1472 | for (d = 0; d <= m; d++) | |||
1473 | { | |||
1474 | fr.x[i][d] += fr.box[m][d]; | |||
1475 | } | |||
1476 | } | |||
1477 | while (fr.x[i][m]-xp[i][m] > hbox[m]) | |||
1478 | { | |||
1479 | for (d = 0; d <= m; d++) | |||
1480 | { | |||
1481 | fr.x[i][d] -= fr.box[m][d]; | |||
1482 | } | |||
1483 | } | |||
1484 | } | |||
1485 | } | |||
1486 | } | |||
1487 | } | |||
1488 | else if (bCluster) | |||
1489 | { | |||
1490 | calc_pbc_cluster(ecenter, ifit, &top, ePBC, fr.x, ind_fit, fr.box); | |||
1491 | } | |||
1492 | ||||
1493 | if (bPFit) | |||
1494 | { | |||
1495 | /* Now modify the coords according to the flags, | |||
1496 | for normal fit, this is only done for output frames */ | |||
1497 | if (bRmPBC) | |||
1498 | { | |||
1499 | gmx_rmpbc_trxfr(gpbc, &fr); | |||
1500 | } | |||
1501 | ||||
1502 | reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL((void*)0), fr.x, w_rls); | |||
1503 | do_fit(natoms, w_rls, xp, fr.x); | |||
1504 | } | |||
1505 | ||||
1506 | /* store this set of coordinates for future use */ | |||
1507 | if (bPFit || bNoJump) | |||
1508 | { | |||
1509 | if (xp == NULL((void*)0)) | |||
1510 | { | |||
1511 | snew(xp, natoms)(xp) = save_calloc("xp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1511, (natoms), sizeof(*(xp))); | |||
1512 | } | |||
1513 | for (i = 0; (i < natoms); i++) | |||
1514 | { | |||
1515 | copy_rvec(fr.x[i], xp[i]); | |||
1516 | rvec_inc(fr.x[i], x_shift); | |||
1517 | } | |||
1518 | } | |||
1519 | ||||
1520 | if (frindex) | |||
1521 | { | |||
1522 | /* see if we have a frame from the frame index group */ | |||
1523 | for (i = 0; i < nrfri && !bDumpFrame; i++) | |||
1524 | { | |||
1525 | bDumpFrame = frame == frindex[i]; | |||
1526 | } | |||
1527 | } | |||
1528 | if (debug && bDumpFrame) | |||
1529 | { | |||
1530 | fprintf(debug, "dumping %d\n", frame); | |||
1531 | } | |||
1532 | ||||
1533 | bWriteFrame = | |||
1534 | ( ( !bTDump && !frindex && frame % skip_nr == 0 ) || bDumpFrame ); | |||
1535 | ||||
1536 | if (bWriteFrame && (bDropUnder || bDropOver)) | |||
1537 | { | |||
1538 | while (dropval[0][drop1] < fr.time && drop1+1 < ndrop) | |||
1539 | { | |||
1540 | drop0 = drop1; | |||
1541 | drop1++; | |||
1542 | } | |||
1543 | if (fabs(dropval[0][drop0] - fr.time) | |||
1544 | < fabs(dropval[0][drop1] - fr.time)) | |||
1545 | { | |||
1546 | dropuse = drop0; | |||
1547 | } | |||
1548 | else | |||
1549 | { | |||
1550 | dropuse = drop1; | |||
1551 | } | |||
1552 | if ((bDropUnder && dropval[1][dropuse] < dropunder) || | |||
1553 | (bDropOver && dropval[1][dropuse] > dropover)) | |||
1554 | { | |||
1555 | bWriteFrame = FALSE0; | |||
1556 | } | |||
1557 | } | |||
1558 | ||||
1559 | if (bWriteFrame) | |||
1560 | { | |||
1561 | ||||
1562 | /* calc new time */ | |||
1563 | if (bTimeStep) | |||
1564 | { | |||
1565 | fr.time = tzero+frame*timestep; | |||
1566 | } | |||
1567 | else | |||
1568 | if (bSetTime) | |||
1569 | { | |||
1570 | fr.time += tshift; | |||
1571 | } | |||
1572 | ||||
1573 | if (bTDump) | |||
1574 | { | |||
1575 | fprintf(stderrstderr, "\nDumping frame at t= %g %s\n", | |||
1576 | output_env_conv_time(oenv, fr.time), output_env_get_time_unit(oenv)); | |||
1577 | } | |||
1578 | ||||
1579 | /* check for writing at each delta_t */ | |||
1580 | bDoIt = (delta_t == 0); | |||
1581 | if (!bDoIt) | |||
1582 | { | |||
1583 | if (!bRound) | |||
1584 | { | |||
1585 | bDoIt = bRmod(fr.time, tzero, delta_t)bRmod_fd(fr.time, tzero, delta_t, 0); | |||
1586 | } | |||
1587 | else | |||
1588 | { | |||
1589 | /* round() is not C89 compatible, so we do this: */ | |||
1590 | bDoIt = bRmod(floor(fr.time+0.5), floor(tzero+0.5),bRmod_fd(floor(fr.time+0.5), floor(tzero+0.5), floor(delta_t+ 0.5), 0) | |||
1591 | floor(delta_t+0.5))bRmod_fd(floor(fr.time+0.5), floor(tzero+0.5), floor(delta_t+ 0.5), 0); | |||
1592 | } | |||
1593 | } | |||
1594 | ||||
1595 | if (bDoIt || bTDump) | |||
1596 | { | |||
1597 | /* print sometimes */ | |||
1598 | if ( ((outframe % SKIP10) == 0) || (outframe < SKIP10) ) | |||
1599 | { | |||
1600 | fprintf(stderrstderr, " -> frame %6d time %8.3f \r", | |||
1601 | outframe, output_env_conv_time(oenv, fr.time)); | |||
1602 | } | |||
1603 | ||||
1604 | if (!bPFit) | |||
1605 | { | |||
1606 | /* Now modify the coords according to the flags, | |||
1607 | for PFit we did this already! */ | |||
1608 | ||||
1609 | if (bRmPBC) | |||
1610 | { | |||
1611 | gmx_rmpbc_trxfr(gpbc, &fr); | |||
1612 | } | |||
1613 | ||||
1614 | if (bReset) | |||
1615 | { | |||
1616 | reset_x_ndim(nfitdim, ifit, ind_fit, natoms, NULL((void*)0), fr.x, w_rls); | |||
1617 | if (bFit) | |||
1618 | { | |||
1619 | do_fit_ndim(nfitdim, natoms, w_rls, xp, fr.x); | |||
1620 | } | |||
1621 | if (!bCenter) | |||
1622 | { | |||
1623 | for (i = 0; i < natoms; i++) | |||
1624 | { | |||
1625 | rvec_inc(fr.x[i], x_shift); | |||
1626 | } | |||
1627 | } | |||
1628 | } | |||
1629 | ||||
1630 | if (bCenter) | |||
1631 | { | |||
1632 | center_x(ecenter, fr.x, fr.box, natoms, ncent, cindex); | |||
1633 | } | |||
1634 | } | |||
1635 | ||||
1636 | if (bPBCcomAtom) | |||
1637 | { | |||
1638 | switch (unitcell_enum) | |||
1639 | { | |||
1640 | case euRect: | |||
1641 | put_atoms_in_box(ePBC, fr.box, natoms, fr.x); | |||
1642 | break; | |||
1643 | case euTric: | |||
1644 | put_atoms_in_triclinic_unitcell(ecenter, fr.box, natoms, fr.x); | |||
1645 | break; | |||
1646 | case euCompact: | |||
1647 | warn = put_atoms_in_compact_unitcell(ePBC, ecenter, fr.box, | |||
1648 | natoms, fr.x); | |||
1649 | if (warn && !bWarnCompact) | |||
1650 | { | |||
1651 | fprintf(stderrstderr, "\n%s\n", warn); | |||
1652 | bWarnCompact = TRUE1; | |||
1653 | } | |||
1654 | break; | |||
1655 | } | |||
1656 | } | |||
1657 | if (bPBCcomRes) | |||
1658 | { | |||
1659 | put_residue_com_in_box(unitcell_enum, ecenter, | |||
1660 | natoms, atoms->atom, ePBC, fr.box, fr.x); | |||
1661 | } | |||
1662 | if (bPBCcomMol) | |||
1663 | { | |||
1664 | put_molecule_com_in_box(unitcell_enum, ecenter, | |||
1665 | &top.mols, | |||
1666 | natoms, atoms->atom, ePBC, fr.box, fr.x); | |||
1667 | } | |||
1668 | /* Copy the input trxframe struct to the output trxframe struct */ | |||
1669 | frout = fr; | |||
1670 | frout.bV = (frout.bV && bVels); | |||
1671 | frout.bF = (frout.bF && bForce); | |||
1672 | frout.natoms = nout; | |||
1673 | if (bNeedPrec && (bSetPrec || !fr.bPrec)) | |||
1674 | { | |||
1675 | frout.bPrec = TRUE1; | |||
1676 | frout.prec = prec; | |||
1677 | } | |||
1678 | if (bCopy) | |||
1679 | { | |||
1680 | frout.x = xmem; | |||
1681 | if (frout.bV) | |||
1682 | { | |||
1683 | frout.v = vmem; | |||
1684 | } | |||
1685 | if (frout.bF) | |||
1686 | { | |||
1687 | frout.f = fmem; | |||
1688 | } | |||
1689 | for (i = 0; i < nout; i++) | |||
1690 | { | |||
1691 | copy_rvec(fr.x[index[i]], frout.x[i]); | |||
1692 | if (frout.bV) | |||
1693 | { | |||
1694 | copy_rvec(fr.v[index[i]], frout.v[i]); | |||
1695 | } | |||
1696 | if (frout.bF) | |||
1697 | { | |||
1698 | copy_rvec(fr.f[index[i]], frout.f[i]); | |||
1699 | } | |||
1700 | } | |||
1701 | } | |||
1702 | ||||
1703 | if (opt2parg_bSet("-shift", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa)) | |||
1704 | { | |||
1705 | for (i = 0; i < nout; i++) | |||
1706 | { | |||
1707 | for (d = 0; d < DIM3; d++) | |||
1708 | { | |||
1709 | frout.x[i][d] += outframe*shift[d]; | |||
1710 | } | |||
1711 | } | |||
1712 | } | |||
1713 | ||||
1714 | if (!bRound) | |||
1715 | { | |||
1716 | bSplitHere = bSplit && bRmod(fr.time, tzero, split_t)bRmod_fd(fr.time, tzero, split_t, 0); | |||
1717 | } | |||
1718 | else | |||
1719 | { | |||
1720 | /* round() is not C89 compatible, so we do this: */ | |||
1721 | bSplitHere = bSplit && bRmod(floor(fr.time+0.5),bRmod_fd(floor(fr.time+0.5), floor(tzero+0.5), floor(split_t+ 0.5), 0) | |||
1722 | floor(tzero+0.5),bRmod_fd(floor(fr.time+0.5), floor(tzero+0.5), floor(split_t+ 0.5), 0) | |||
1723 | floor(split_t+0.5))bRmod_fd(floor(fr.time+0.5), floor(tzero+0.5), floor(split_t+ 0.5), 0); | |||
1724 | } | |||
1725 | if (bSeparate || bSplitHere) | |||
1726 | { | |||
1727 | mk_filenm(outf_base, ftp2ext(ftp), nzero, file_nr, out_file2); | |||
1728 | } | |||
1729 | ||||
1730 | switch (ftp) | |||
1731 | { | |||
1732 | case efTNG: | |||
1733 | write_tng_frame(trxout, &frout); | |||
1734 | // TODO when trjconv behaves better: work how to read and write lambda | |||
1735 | break; | |||
1736 | case efTRJ: | |||
1737 | case efTRR: | |||
1738 | case efXTC: | |||
1739 | if (bSplitHere) | |||
1740 | { | |||
1741 | if (trxout) | |||
1742 | { | |||
1743 | close_trx(trxout); | |||
1744 | } | |||
1745 | trxout = open_trx(out_file2, filemode); | |||
1746 | } | |||
1747 | if (bSubTraj) | |||
1748 | { | |||
1749 | if (my_clust != -1) | |||
1750 | { | |||
1751 | char buf[STRLEN4096]; | |||
1752 | if (clust_status_id[my_clust] == -1) | |||
1753 | { | |||
1754 | sprintf(buf, "%s.%s", clust->grpname[my_clust], ftp2ext(ftp)); | |||
1755 | clust_status[my_clust] = open_trx(buf, "w"); | |||
1756 | clust_status_id[my_clust] = 1; | |||
1757 | ntrxopen++; | |||
1758 | } | |||
1759 | else if (clust_status_id[my_clust] == -2) | |||
1760 | { | |||
1761 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1761, "File %s.xtc should still be open (%d open .xtc files)\n" "in order to write frame %d. my_clust = %d", | |||
1762 | clust->grpname[my_clust], ntrxopen, frame, | |||
1763 | my_clust); | |||
1764 | } | |||
1765 | write_trxframe(clust_status[my_clust], &frout, gc); | |||
1766 | nfwritten[my_clust]++; | |||
1767 | if (nfwritten[my_clust] == | |||
1768 | (clust->clust->index[my_clust+1]- | |||
1769 | clust->clust->index[my_clust])) | |||
1770 | { | |||
1771 | close_trx(clust_status[my_clust]); | |||
1772 | clust_status[my_clust] = NULL((void*)0); | |||
1773 | clust_status_id[my_clust] = -2; | |||
1774 | ntrxopen--; | |||
1775 | if (ntrxopen < 0) | |||
1776 | { | |||
1777 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1777, "Less than zero open .xtc files!"); | |||
1778 | } | |||
1779 | } | |||
1780 | } | |||
1781 | } | |||
1782 | else | |||
1783 | { | |||
1784 | write_trxframe(trxout, &frout, gc); | |||
1785 | } | |||
1786 | break; | |||
1787 | case efGRO: | |||
1788 | case efG96: | |||
1789 | case efPDB: | |||
1790 | sprintf(title, "Generated by trjconv : %s t= %9.5f", | |||
1791 | top_title, fr.time); | |||
1792 | if (bSeparate || bSplitHere) | |||
1793 | { | |||
1794 | out = gmx_ffopen(out_file2, "w"); | |||
1795 | } | |||
1796 | switch (ftp) | |||
1797 | { | |||
1798 | case efGRO: | |||
1799 | write_hconf_p(out, title, &useatoms, prec2ndec(frout.prec), | |||
1800 | frout.x, frout.bV ? frout.v : NULL((void*)0), frout.box); | |||
1801 | break; | |||
1802 | case efPDB: | |||
1803 | fprintf(out, "REMARK GENERATED BY TRJCONV\n"); | |||
1804 | sprintf(title, "%s t= %9.5f", top_title, frout.time); | |||
1805 | /* if reading from pdb, we want to keep the original | |||
1806 | model numbering else we write the output frame | |||
1807 | number plus one, because model 0 is not allowed in pdb */ | |||
1808 | if (ftpin == efPDB && fr.bStep && fr.step > model_nr) | |||
1809 | { | |||
1810 | model_nr = fr.step; | |||
1811 | } | |||
1812 | else | |||
1813 | { | |||
1814 | model_nr++; | |||
1815 | } | |||
1816 | write_pdbfile(out, title, &useatoms, frout.x, | |||
1817 | frout.ePBC, frout.box, ' ', model_nr, gc, TRUE1); | |||
1818 | break; | |||
1819 | case efG96: | |||
1820 | frout.title = title; | |||
1821 | if (bSeparate || bTDump) | |||
1822 | { | |||
1823 | frout.bTitle = TRUE1; | |||
1824 | if (bTPS) | |||
1825 | { | |||
1826 | frout.bAtoms = TRUE1; | |||
1827 | } | |||
1828 | frout.atoms = &useatoms; | |||
1829 | frout.bStep = FALSE0; | |||
1830 | frout.bTime = FALSE0; | |||
1831 | } | |||
1832 | else | |||
1833 | { | |||
1834 | frout.bTitle = (outframe == 0); | |||
1835 | frout.bAtoms = FALSE0; | |||
1836 | frout.bStep = TRUE1; | |||
1837 | frout.bTime = TRUE1; | |||
1838 | } | |||
1839 | write_g96_conf(out, &frout, -1, NULL((void*)0)); | |||
1840 | } | |||
1841 | if (bSeparate) | |||
1842 | { | |||
1843 | gmx_ffclose(out); | |||
1844 | out = NULL((void*)0); | |||
1845 | } | |||
1846 | break; | |||
1847 | default: | |||
1848 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1848, "DHE, ftp=%d\n", ftp); | |||
1849 | } | |||
1850 | if (bSeparate || bSplitHere) | |||
1851 | { | |||
1852 | file_nr++; | |||
1853 | } | |||
1854 | ||||
1855 | /* execute command */ | |||
1856 | if (bExec) | |||
1857 | { | |||
1858 | char c[255]; | |||
1859 | sprintf(c, "%s %d", exec_command, file_nr-1); | |||
1860 | /*fprintf(stderr,"Executing '%s'\n",c);*/ | |||
1861 | #ifdef GMX_NO_SYSTEM | |||
1862 | printf("Warning-- No calls to system(3) supported on this platform."); | |||
1863 | printf("Warning-- Skipping execution of 'system(\"%s\")'.", c); | |||
1864 | #else | |||
1865 | if (0 != system(c)) | |||
1866 | { | |||
1867 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1867, "Error executing command: %s", c); | |||
1868 | } | |||
1869 | #endif | |||
1870 | } | |||
1871 | outframe++; | |||
1872 | } | |||
1873 | } | |||
1874 | frame++; | |||
1875 | bHaveNextFrame = read_next_frame(oenv, trxin, &fr); | |||
1876 | } | |||
1877 | while (!(bTDump && bDumpFrame) && bHaveNextFrame); | |||
1878 | } | |||
1879 | ||||
1880 | if (!bHaveFirstFrame || (bTDump && !bDumpFrame)) | |||
1881 | { | |||
1882 | fprintf(stderrstderr, "\nWARNING no output, " | |||
1883 | "last frame read at t=%g\n", fr.time); | |||
1884 | } | |||
1885 | fprintf(stderrstderr, "\n"); | |||
1886 | ||||
1887 | close_trj(trxin); | |||
1888 | sfree(outf_base)save_free("outf_base", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1888, (outf_base)); | |||
1889 | ||||
1890 | if (bRmPBC) | |||
1891 | { | |||
1892 | gmx_rmpbc_done(gpbc); | |||
1893 | } | |||
1894 | ||||
1895 | if (trxout) | |||
1896 | { | |||
1897 | close_trx(trxout); | |||
1898 | } | |||
1899 | else if (out != NULL((void*)0)) | |||
1900 | { | |||
1901 | gmx_ffclose(out); | |||
1902 | } | |||
1903 | if (bSubTraj) | |||
1904 | { | |||
1905 | for (i = 0; (i < clust->clust->nr); i++) | |||
1906 | { | |||
1907 | if (clust_status_id[i] >= 0) | |||
1908 | { | |||
1909 | close_trx(clust_status[i]); | |||
1910 | } | |||
1911 | } | |||
1912 | } | |||
1913 | } | |||
1914 | ||||
1915 | sfree(mtop)save_free("mtop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjconv.c" , 1915, (mtop)); | |||
1916 | ||||
1917 | do_view(oenv, out_file, NULL((void*)0)); | |||
1918 | ||||
1919 | return 0; | |||
1920 | } |