File: | gromacs/gmxana/gmx_trjconv.c |
Location: | line 1000, column 17 |
Description: | Value stored to 'bSetUR' is never read |
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; |
Value stored to 'bSetUR' is never read | |
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 | } |