2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, 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.
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.
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.
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.
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.
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.
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/commandline/viewit.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxana/gmx_ana.h"
52 #include "gromacs/linearalgebra/nrjac.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/mdtypes/md_enums.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/pbcutil/rmpbc.h"
59 #include "gromacs/topology/index.h"
60 #include "gromacs/topology/topology.h"
61 #include "gromacs/trajectory/trajectoryframe.h"
62 #include "gromacs/utility/arraysize.h"
63 #include "gromacs/utility/cstringutil.h"
64 #include "gromacs/utility/fatalerror.h"
65 #include "gromacs/utility/futil.h"
66 #include "gromacs/utility/smalloc.h"
68 static void low_print_data(FILE *fp
, real time
, rvec x
[], int n
, int *index
,
69 gmx_bool bDim
[], const char *sffmt
)
73 fprintf(fp
, " %g", time
);
74 for (i
= 0; i
< n
; i
++)
84 for (d
= 0; d
< DIM
; d
++)
88 fprintf(fp
, sffmt
, x
[ii
][d
]);
93 fprintf(fp
, sffmt
, norm(x
[ii
]));
99 static void average_data(rvec x
[], rvec xav
[], real
*mass
,
100 int ngrps
, int isize
[], int **index
)
105 double sum
[DIM
], mtot
;
107 for (g
= 0; g
< ngrps
; g
++)
112 for (i
= 0; i
< isize
[g
]; i
++)
118 svmul(m
, x
[ind
], tmp
);
119 for (d
= 0; d
< DIM
; d
++)
127 for (d
= 0; d
< DIM
; d
++)
135 for (d
= 0; d
< DIM
; d
++)
137 xav
[g
][d
] = sum
[d
]/mtot
;
142 /* mass=NULL, so these are forces: sum only (do not average) */
143 for (d
= 0; d
< DIM
; d
++)
151 static void print_data(FILE *fp
, real time
, rvec x
[], real
*mass
, gmx_bool bCom
,
152 int ngrps
, int isize
[], int **index
, gmx_bool bDim
[],
155 static rvec
*xav
= nullptr;
163 average_data(x
, xav
, mass
, ngrps
, isize
, index
);
164 low_print_data(fp
, time
, xav
, ngrps
, nullptr, bDim
, sffmt
);
168 low_print_data(fp
, time
, x
, isize
[0], index
[0], bDim
, sffmt
);
172 static void write_trx_x(t_trxstatus
*status
, const t_trxframe
*fr
, real
*mass
, gmx_bool bCom
,
173 int ngrps
, int isize
[], int **index
)
175 static rvec
*xav
= nullptr;
176 static t_atoms
*atoms
= nullptr;
187 snew(atoms
->atom
, ngrps
);
189 /* Note that atom and residue names will be the ones
190 * of the first atom in each group.
192 for (i
= 0; i
< ngrps
; i
++)
194 atoms
->atom
[i
] = fr
->atoms
->atom
[index
[i
][0]];
195 atoms
->atomname
[i
] = fr
->atoms
->atomname
[index
[i
][0]];
198 average_data(fr
->x
, xav
, mass
, ngrps
, isize
, index
);
200 fr_av
.natoms
= ngrps
;
203 write_trxframe(status
, &fr_av
, nullptr);
207 write_trxframe_indexed(status
, fr
, isize
[0], index
[0], nullptr);
211 static void make_legend(FILE *fp
, int ngrps
, int isize
, int index
[],
212 char **name
, gmx_bool bCom
, gmx_bool bMol
, gmx_bool bDim
[],
213 const gmx_output_env_t
*oenv
)
216 const char *dimtxt
[] = { " X", " Y", " Z", "" };
230 for (i
= 0; i
< n
; i
++)
232 for (d
= 0; d
<= DIM
; d
++)
236 snew(leg
[j
], STRLEN
);
239 sprintf(leg
[j
], "mol %d%s", index
[i
]+1, dimtxt
[d
]);
243 sprintf(leg
[j
], "%s%s", name
[i
], dimtxt
[d
]);
247 sprintf(leg
[j
], "atom %d%s", index
[i
]+1, dimtxt
[d
]);
253 xvgr_legend(fp
, j
, (const char**)leg
, oenv
);
255 for (i
= 0; i
< j
; i
++)
262 static real
ekrot(rvec x
[], rvec v
[], real mass
[], int isize
, int index
[])
264 static real
**TCM
= nullptr, **L
;
265 double tm
, m0
, lxx
, lxy
, lxz
, lyy
, lyz
, lzz
, ekrot
;
274 for (i
= 0; i
< DIM
; i
++)
279 for (i
= 0; i
< DIM
; i
++)
289 for (i
= 0; i
< isize
; i
++)
294 cprod(x
[j
], v
[j
], a0
);
295 for (m
= 0; (m
< DIM
); m
++)
297 xcm
[m
] += m0
*x
[j
][m
]; /* c.o.m. position */
298 vcm
[m
] += m0
*v
[j
][m
]; /* c.o.m. velocity */
299 acm
[m
] += m0
*a0
[m
]; /* rotational velocity around c.o.m. */
302 dcprod(xcm
, vcm
, b0
);
303 for (m
= 0; (m
< DIM
); m
++)
310 lxx
= lxy
= lxz
= lyy
= lyz
= lzz
= 0.0;
311 for (i
= 0; i
< isize
; i
++)
315 for (m
= 0; m
< DIM
; m
++)
317 dx
[m
] = x
[j
][m
] - xcm
[m
];
319 lxx
+= dx
[XX
]*dx
[XX
]*m0
;
320 lxy
+= dx
[XX
]*dx
[YY
]*m0
;
321 lxz
+= dx
[XX
]*dx
[ZZ
]*m0
;
322 lyy
+= dx
[YY
]*dx
[YY
]*m0
;
323 lyz
+= dx
[YY
]*dx
[ZZ
]*m0
;
324 lzz
+= dx
[ZZ
]*dx
[ZZ
]*m0
;
327 L
[XX
][XX
] = lyy
+ lzz
;
331 L
[YY
][YY
] = lxx
+ lzz
;
335 L
[ZZ
][ZZ
] = lxx
+ lyy
;
337 m_inv_gen(L
, DIM
, TCM
);
339 /* Compute omega (hoeksnelheid) */
342 for (m
= 0; m
< DIM
; m
++)
344 for (n
= 0; n
< DIM
; n
++)
346 ocm
[m
] += TCM
[m
][n
]*acm
[n
];
348 ekrot
+= 0.5*ocm
[m
]*acm
[m
];
354 static real
ektrans(rvec v
[], real mass
[], int isize
, int index
[])
362 for (i
= 0; i
< isize
; i
++)
365 for (d
= 0; d
< DIM
; d
++)
367 mvcom
[d
] += mass
[j
]*v
[j
][d
];
372 return dnorm2(mvcom
)/(mtot
*2);
375 static real
temp(rvec v
[], real mass
[], int isize
, int index
[])
381 for (i
= 0; i
< isize
; i
++)
384 ekin2
+= mass
[j
]*norm2(v
[j
]);
387 return ekin2
/(3*isize
*BOLTZ
);
390 static void remove_jump(matrix box
, int natoms
, rvec xp
[], rvec x
[])
395 for (d
= 0; d
< DIM
; d
++)
397 hbox
[d
] = 0.5*box
[d
][d
];
399 for (i
= 0; i
< natoms
; i
++)
401 for (m
= DIM
-1; m
>= 0; m
--)
403 while (x
[i
][m
]-xp
[i
][m
] <= -hbox
[m
])
405 for (d
= 0; d
<= m
; d
++)
407 x
[i
][d
] += box
[m
][d
];
410 while (x
[i
][m
]-xp
[i
][m
] > hbox
[m
])
412 for (d
= 0; d
<= m
; d
++)
414 x
[i
][d
] -= box
[m
][d
];
421 static void write_pdb_bfac(const char *fname
, const char *xname
,
422 const char *title
, t_atoms
*atoms
, int ePBC
, matrix box
,
423 int isize
, int *index
, int nfr_x
, rvec
*x
,
424 int nfr_v
, rvec
*sum
,
425 gmx_bool bDim
[], real scale_factor
,
426 const gmx_output_env_t
*oenv
)
429 real max
, len2
, scale
;
433 if ((nfr_x
== 0) || (nfr_v
== 0))
435 fprintf(stderr
, "No frames found for %s, will not write %s\n",
440 fprintf(stderr
, "Used %d frames for %s\n", nfr_x
, "coordinates");
441 fprintf(stderr
, "Used %d frames for %s\n", nfr_v
, title
);
446 for (i
= 0; i
< DIM
; i
++)
460 for (i
= 0; i
< isize
; i
++)
462 svmul(scale
, sum
[index
[i
]], sum
[index
[i
]]);
465 fp
= xvgropen(xname
, title
, "Atom", "Spatial component", oenv
);
466 for (i
= 0; i
< isize
; i
++)
468 fprintf(fp
, "%-5d %10.3f %10.3f %10.3f\n", 1+i
,
469 sum
[index
[i
]][XX
], sum
[index
[i
]][YY
], sum
[index
[i
]][ZZ
]);
474 for (i
= 0; i
< isize
; i
++)
477 for (m
= 0; m
< DIM
; m
++)
479 if (bDim
[m
] || bDim
[DIM
])
481 len2
+= gmx::square(sum
[index
[i
]][m
]);
490 if (scale_factor
!= 0)
492 scale
= scale_factor
;
502 scale
= 10.0/std::sqrt(max
);
506 printf("Maximum %s is %g on atom %d %s, res. %s %d\n",
507 title
, std::sqrt(max
), maxi
+1, *(atoms
->atomname
[maxi
]),
508 *(atoms
->resinfo
[atoms
->atom
[maxi
].resind
].name
),
509 atoms
->resinfo
[atoms
->atom
[maxi
].resind
].nr
);
511 if (atoms
->pdbinfo
== nullptr)
513 snew(atoms
->pdbinfo
, atoms
->nr
);
515 atoms
->havePdbInfo
= TRUE
;
519 for (i
= 0; i
< isize
; i
++)
522 for (m
= 0; m
< DIM
; m
++)
524 if (bDim
[m
] || bDim
[DIM
])
526 len2
+= gmx::square(sum
[index
[i
]][m
]);
529 atoms
->pdbinfo
[index
[i
]].bfac
= std::sqrt(len2
)*scale
;
534 for (i
= 0; i
< isize
; i
++)
536 atoms
->pdbinfo
[index
[i
]].bfac
= sum
[index
[i
]][onedim
]*scale
;
539 write_sto_conf_indexed(fname
, title
, atoms
, x
, nullptr, ePBC
, box
, isize
, index
);
543 static void update_histo(int gnx
, int index
[], rvec v
[],
544 int *nhisto
, int **histo
, real binwidth
)
549 if (*histo
== nullptr)
552 for (i
= 0; (i
< gnx
); i
++)
554 vn
= norm(v
[index
[i
]]);
555 vnmax
= std::max(vn
, vnmax
);
558 *nhisto
= static_cast<int>(1+(vnmax
/binwidth
));
559 snew(*histo
, *nhisto
);
561 for (i
= 0; (i
< gnx
); i
++)
563 vn
= norm(v
[index
[i
]]);
564 in
= static_cast<int>(vn
/binwidth
);
568 fprintf(stderr
, "Extending histogram from %d to %d\n", *nhisto
, nnn
);
571 for (m
= *nhisto
; (m
< nnn
); m
++)
581 static void print_histo(const char *fn
, int nhisto
, int histo
[], real binwidth
,
582 const gmx_output_env_t
*oenv
)
587 fp
= xvgropen(fn
, "Velocity distribution", "V (nm/ps)", "arbitrary units",
589 for (i
= 0; (i
< nhisto
); i
++)
591 fprintf(fp
, "%10.3e %10d\n", i
*binwidth
, histo
[i
]);
596 int gmx_traj(int argc
, char *argv
[])
598 const char *desc
[] = {
599 "[THISMODULE] plots coordinates, velocities, forces and/or the box.",
600 "With [TT]-com[tt] the coordinates, velocities and forces are",
601 "calculated for the center of mass of each group.",
602 "When [TT]-mol[tt] is set, the numbers in the index file are",
603 "interpreted as molecule numbers and the same procedure as with",
604 "[TT]-com[tt] is used for each molecule.[PAR]",
605 "Option [TT]-ot[tt] plots the temperature of each group,",
606 "provided velocities are present in the trajectory file.",
607 "No corrections are made for constrained degrees of freedom!",
608 "This implies [TT]-com[tt].[PAR]",
609 "Options [TT]-ekt[tt] and [TT]-ekr[tt] plot the translational and",
610 "rotational kinetic energy of each group,",
611 "provided velocities are present in the trajectory file.",
612 "This implies [TT]-com[tt].[PAR]",
613 "Options [TT]-cv[tt] and [TT]-cf[tt] write the average velocities",
614 "and average forces as temperature factors to a [REF].pdb[ref] file with",
615 "the average coordinates or the coordinates at [TT]-ctime[tt].",
616 "The temperature factors are scaled such that the maximum is 10.",
617 "The scaling can be changed with the option [TT]-scale[tt].",
618 "To get the velocities or forces of one",
619 "frame set both [TT]-b[tt] and [TT]-e[tt] to the time of",
620 "desired frame. When averaging over frames you might need to use",
621 "the [TT]-nojump[tt] option to obtain the correct average coordinates.",
622 "If you select either of these option the average force and velocity",
623 "for each atom are written to an [REF].xvg[ref] file as well",
624 "(specified with [TT]-av[tt] or [TT]-af[tt]).[PAR]",
625 "Option [TT]-vd[tt] computes a velocity distribution, i.e. the",
626 "norm of the vector is plotted. In addition in the same graph",
627 "the kinetic energy distribution is given.",
629 "See [gmx-trajectory] for plotting similar data for selections."
631 static gmx_bool bMol
= FALSE
, bCom
= FALSE
, bPBC
= TRUE
, bNoJump
= FALSE
;
632 static gmx_bool bX
= TRUE
, bY
= TRUE
, bZ
= TRUE
, bNorm
= FALSE
, bFP
= FALSE
;
633 static int ngroups
= 1;
634 static real ctime
= -1, scale
= 0, binwidth
= 1;
636 { "-com", FALSE
, etBOOL
, {&bCom
},
637 "Plot data for the com of each group" },
638 { "-pbc", FALSE
, etBOOL
, {&bPBC
},
639 "Make molecules whole for COM" },
640 { "-mol", FALSE
, etBOOL
, {&bMol
},
641 "Index contains molecule numbers iso atom numbers" },
642 { "-nojump", FALSE
, etBOOL
, {&bNoJump
},
643 "Remove jumps of atoms across the box" },
644 { "-x", FALSE
, etBOOL
, {&bX
},
645 "Plot X-component" },
646 { "-y", FALSE
, etBOOL
, {&bY
},
647 "Plot Y-component" },
648 { "-z", FALSE
, etBOOL
, {&bZ
},
649 "Plot Z-component" },
650 { "-ng", FALSE
, etINT
, {&ngroups
},
651 "Number of groups to consider" },
652 { "-len", FALSE
, etBOOL
, {&bNorm
},
653 "Plot vector length" },
654 { "-fp", FALSE
, etBOOL
, {&bFP
},
655 "Full precision output" },
656 { "-bin", FALSE
, etREAL
, {&binwidth
},
657 "Binwidth for velocity histogram (nm/ps)" },
658 { "-ctime", FALSE
, etREAL
, {&ctime
},
659 "Use frame at this time for x in [TT]-cv[tt] and [TT]-cf[tt] instead of the average x" },
660 { "-scale", FALSE
, etREAL
, {&scale
},
661 "Scale factor for [REF].pdb[ref] output, 0 is autoscale" }
663 FILE *outx
= nullptr, *outv
= nullptr, *outf
= nullptr, *outb
= nullptr, *outt
= nullptr;
664 FILE *outekt
= nullptr, *outekr
= nullptr;
670 int flags
, nvhisto
= 0, *vhisto
= nullptr;
671 rvec
*xtop
, *xp
= nullptr;
672 rvec
*sumx
= nullptr, *sumv
= nullptr, *sumf
= nullptr;
675 t_trxstatus
*status_out
= nullptr;
676 gmx_rmpbc_t gpbc
= nullptr;
678 int nr_xfr
, nr_vfr
, nr_ffr
;
681 int **index0
, **index
;
684 gmx_bool bTop
, bOX
, bOXT
, bOV
, bOF
, bOB
, bOT
, bEKT
, bEKR
, bCV
, bCF
;
685 gmx_bool bDim
[4], bDum
[4], bVD
;
686 char sffmt
[STRLEN
], sffmt6
[STRLEN
];
687 const char *box_leg
[6] = { "XX", "YY", "ZZ", "YX", "ZX", "ZY" };
688 gmx_output_env_t
*oenv
;
691 { efTRX
, "-f", nullptr, ffREAD
},
692 { efTPS
, nullptr, nullptr, ffREAD
},
693 { efNDX
, nullptr, nullptr, ffOPTRD
},
694 { efXVG
, "-ox", "coord", ffOPTWR
},
695 { efTRX
, "-oxt", "coord", ffOPTWR
},
696 { efXVG
, "-ov", "veloc", ffOPTWR
},
697 { efXVG
, "-of", "force", ffOPTWR
},
698 { efXVG
, "-ob", "box", ffOPTWR
},
699 { efXVG
, "-ot", "temp", ffOPTWR
},
700 { efXVG
, "-ekt", "ektrans", ffOPTWR
},
701 { efXVG
, "-ekr", "ekrot", ffOPTWR
},
702 { efXVG
, "-vd", "veldist", ffOPTWR
},
703 { efPDB
, "-cv", "veloc", ffOPTWR
},
704 { efPDB
, "-cf", "force", ffOPTWR
},
705 { efXVG
, "-av", "all_veloc", ffOPTWR
},
706 { efXVG
, "-af", "all_force", ffOPTWR
}
708 #define NFILE asize(fnm)
710 if (!parse_common_args(&argc
, argv
,
711 PCA_CAN_TIME
| PCA_TIME_UNIT
| PCA_CAN_VIEW
,
712 NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
, 0, nullptr, &oenv
))
719 fprintf(stderr
, "Interpreting indexfile entries as molecules.\n"
720 "Using center of mass.\n");
723 bOX
= opt2bSet("-ox", NFILE
, fnm
);
724 bOXT
= opt2bSet("-oxt", NFILE
, fnm
);
725 bOV
= opt2bSet("-ov", NFILE
, fnm
);
726 bOF
= opt2bSet("-of", NFILE
, fnm
);
727 bOB
= opt2bSet("-ob", NFILE
, fnm
);
728 bOT
= opt2bSet("-ot", NFILE
, fnm
);
729 bEKT
= opt2bSet("-ekt", NFILE
, fnm
);
730 bEKR
= opt2bSet("-ekr", NFILE
, fnm
);
731 bCV
= opt2bSet("-cv", NFILE
, fnm
) || opt2bSet("-av", NFILE
, fnm
);
732 bCF
= opt2bSet("-cf", NFILE
, fnm
) || opt2bSet("-af", NFILE
, fnm
);
733 bVD
= opt2bSet("-vd", NFILE
, fnm
) || opt2parg_bSet("-bin", asize(pa
), pa
);
734 if (bMol
|| bOT
|| bEKT
|| bEKR
)
746 sprintf(sffmt
, "\t%s", gmx_real_fullprecision_pfmt
);
750 sprintf(sffmt
, "\t%%g");
752 sprintf(sffmt6
, "%s%s%s%s%s%s", sffmt
, sffmt
, sffmt
, sffmt
, sffmt
, sffmt
);
754 bTop
= read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), &top
, &ePBC
,
755 &xtop
, nullptr, topbox
,
756 bCom
&& (bOX
|| bOXT
|| bOV
|| bOT
|| bEKT
|| bEKR
));
758 if ((bMol
|| bCV
|| bCF
) && !bTop
)
760 gmx_fatal(FARGS
, "Need a run input file for option -mol, -cv or -cf");
765 indexfn
= ftp2fn(efNDX
, NFILE
, fnm
);
769 indexfn
= ftp2fn_null(efNDX
, NFILE
, fnm
);
772 if (!(bCom
&& !bMol
))
776 snew(grpname
, ngroups
);
777 snew(isize0
, ngroups
);
778 snew(index0
, ngroups
);
779 get_index(&(top
.atoms
), indexfn
, ngroups
, isize0
, index0
, grpname
);
786 snew(isize
, ngroups
);
787 snew(index
, ngroups
);
788 for (i
= 0; i
< ngroups
; i
++)
790 if (index0
[0][i
] < 0 || index0
[0][i
] >= mols
->nr
)
792 gmx_fatal(FARGS
, "Molecule index (%d) is out of range (%d-%d)",
793 index0
[0][i
]+1, 1, mols
->nr
);
795 isize
[i
] = atndx
[index0
[0][i
]+1] - atndx
[index0
[0][i
]];
796 snew(index
[i
], isize
[i
]);
797 for (j
= 0; j
< isize
[i
]; j
++)
799 index
[i
][j
] = atndx
[index0
[0][i
]] + j
;
810 snew(mass
, top
.atoms
.nr
);
811 for (i
= 0; i
< top
.atoms
.nr
; i
++)
813 mass
[i
] = top
.atoms
.atom
[i
].m
;
822 std::string
label(output_env_get_xvgr_tlabel(oenv
));
825 flags
= flags
| TRX_READ_X
;
826 outx
= xvgropen(opt2fn("-ox", NFILE
, fnm
),
827 bCom
? "Center of mass" : "Coordinate",
828 label
, "Coordinate (nm)", oenv
);
829 make_legend(outx
, ngroups
, isize0
[0], index0
[0], grpname
, bCom
, bMol
, bDim
, oenv
);
833 flags
= flags
| TRX_READ_X
;
834 status_out
= open_trx(opt2fn("-oxt", NFILE
, fnm
), "w");
838 flags
= flags
| TRX_READ_V
;
839 outv
= xvgropen(opt2fn("-ov", NFILE
, fnm
),
840 bCom
? "Center of mass velocity" : "Velocity",
841 label
, "Velocity (nm/ps)", oenv
);
842 make_legend(outv
, ngroups
, isize0
[0], index0
[0], grpname
, bCom
, bMol
, bDim
, oenv
);
846 flags
= flags
| TRX_READ_F
;
847 outf
= xvgropen(opt2fn("-of", NFILE
, fnm
), "Force",
848 label
, "Force (kJ mol\\S-1\\N nm\\S-1\\N)",
850 make_legend(outf
, ngroups
, isize0
[0], index0
[0], grpname
, bCom
, bMol
, bDim
, oenv
);
854 outb
= xvgropen(opt2fn("-ob", NFILE
, fnm
), "Box vector elements",
855 label
, "(nm)", oenv
);
857 xvgr_legend(outb
, 6, box_leg
, oenv
);
865 flags
= flags
| TRX_READ_V
;
866 outt
= xvgropen(opt2fn("-ot", NFILE
, fnm
), "Temperature",
868 make_legend(outt
, ngroups
, isize
[0], index
[0], grpname
, bCom
, bMol
, bDum
, oenv
);
876 flags
= flags
| TRX_READ_V
;
877 outekt
= xvgropen(opt2fn("-ekt", NFILE
, fnm
), "Center of mass translation",
878 label
, "Energy (kJ mol\\S-1\\N)", oenv
);
879 make_legend(outekt
, ngroups
, isize
[0], index
[0], grpname
, bCom
, bMol
, bDum
, oenv
);
887 flags
= flags
| TRX_READ_X
| TRX_READ_V
;
888 outekr
= xvgropen(opt2fn("-ekr", NFILE
, fnm
), "Center of mass rotation",
889 label
, "Energy (kJ mol\\S-1\\N)", oenv
);
890 make_legend(outekr
, ngroups
, isize
[0], index
[0], grpname
, bCom
, bMol
, bDum
, oenv
);
894 flags
= flags
| TRX_READ_V
;
898 flags
= flags
| TRX_READ_X
| TRX_READ_V
;
902 flags
= flags
| TRX_READ_X
| TRX_READ_F
;
904 if ((flags
== 0) && !bOB
)
906 fprintf(stderr
, "Please select one or more output file options\n");
910 read_first_frame(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &fr
, flags
);
913 if ((bOV
|| bOF
) && fn2ftp(ftp2fn(efTRX
, NFILE
, fnm
)) == efXTC
)
915 gmx_fatal(FARGS
, "Cannot extract velocities or forces since your input XTC file does not contain them.");
920 snew(sumx
, fr
.natoms
);
924 snew(sumv
, fr
.natoms
);
928 snew(sumf
, fr
.natoms
);
936 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, fr
.natoms
);
941 time
= output_env_conv_time(oenv
, fr
.time
);
943 if (fr
.bX
&& bNoJump
&& fr
.bBox
)
947 remove_jump(fr
.box
, fr
.natoms
, xp
, fr
.x
);
953 for (i
= 0; i
< fr
.natoms
; i
++)
955 copy_rvec(fr
.x
[i
], xp
[i
]);
959 if (fr
.bX
&& bCom
&& bPBC
)
961 gmx_rmpbc_trxfr(gpbc
, &fr
);
966 update_histo(isize
[0], index
[0], fr
.v
, &nvhisto
, &vhisto
, binwidth
);
971 print_data(outx
, time
, fr
.x
, mass
, bCom
, ngroups
, isize
, index
, bDim
, sffmt
);
975 t_trxframe frout
= fr
;
978 frout
.atoms
= &top
.atoms
;
983 write_trx_x(status_out
, &frout
, mass
, bCom
, ngroups
, isize
, index
);
987 print_data(outv
, time
, fr
.v
, mass
, bCom
, ngroups
, isize
, index
, bDim
, sffmt
);
991 print_data(outf
, time
, fr
.f
, nullptr, bCom
, ngroups
, isize
, index
, bDim
, sffmt
);
995 fprintf(outb
, "\t%g", fr
.time
);
996 fprintf(outb
, sffmt6
,
997 fr
.box
[XX
][XX
], fr
.box
[YY
][YY
], fr
.box
[ZZ
][ZZ
],
998 fr
.box
[YY
][XX
], fr
.box
[ZZ
][XX
], fr
.box
[ZZ
][YY
]);
1003 fprintf(outt
, " %g", time
);
1004 for (i
= 0; i
< ngroups
; i
++)
1006 fprintf(outt
, sffmt
, temp(fr
.v
, mass
, isize
[i
], index
[i
]));
1008 fprintf(outt
, "\n");
1012 fprintf(outekt
, " %g", time
);
1013 for (i
= 0; i
< ngroups
; i
++)
1015 fprintf(outekt
, sffmt
, ektrans(fr
.v
, mass
, isize
[i
], index
[i
]));
1017 fprintf(outekt
, "\n");
1019 if (bEKR
&& fr
.bX
&& fr
.bV
)
1021 fprintf(outekr
, " %g", time
);
1022 for (i
= 0; i
< ngroups
; i
++)
1024 fprintf(outekr
, sffmt
, ekrot(fr
.x
, fr
.v
, mass
, isize
[i
], index
[i
]));
1026 fprintf(outekr
, "\n");
1028 if ((bCV
|| bCF
) && fr
.bX
&&
1029 (ctime
< 0 || (fr
.time
>= ctime
*0.999999 &&
1030 fr
.time
<= ctime
*1.000001)))
1032 for (i
= 0; i
< fr
.natoms
; i
++)
1034 rvec_inc(sumx
[i
], fr
.x
[i
]);
1040 for (i
= 0; i
< fr
.natoms
; i
++)
1042 rvec_inc(sumv
[i
], fr
.v
[i
]);
1048 for (i
= 0; i
< fr
.natoms
; i
++)
1050 rvec_inc(sumf
[i
], fr
.f
[i
]);
1056 while (read_next_frame(oenv
, status
, &fr
));
1058 if (gpbc
!= nullptr)
1060 gmx_rmpbc_done(gpbc
);
1063 /* clean up a bit */
1072 close_trx(status_out
);
1101 print_histo(opt2fn("-vd", NFILE
, fnm
), nvhisto
, vhisto
, binwidth
, oenv
);
1108 if (ePBC
!= epbcNONE
&& !bNoJump
)
1110 fprintf(stderr
, "\nWARNING: More than one frame was used for option -cv or -cf\n"
1111 "If atoms jump across the box you should use the -nojump or -ctime option\n\n");
1113 for (i
= 0; i
< isize
[0]; i
++)
1115 svmul(1.0/nr_xfr
, sumx
[index
[0][i
]], sumx
[index
[0][i
]]);
1118 else if (nr_xfr
== 0)
1120 fprintf(stderr
, "\nWARNING: No coordinate frames found for option -cv or -cf\n\n");
1125 write_pdb_bfac(opt2fn("-cv", NFILE
, fnm
),
1126 opt2fn("-av", NFILE
, fnm
), "average velocity", &(top
.atoms
),
1127 ePBC
, topbox
, isize
[0], index
[0], nr_xfr
, sumx
,
1128 nr_vfr
, sumv
, bDim
, scale
, oenv
);
1132 write_pdb_bfac(opt2fn("-cf", NFILE
, fnm
),
1133 opt2fn("-af", NFILE
, fnm
), "average force", &(top
.atoms
),
1134 ePBC
, topbox
, isize
[0], index
[0], nr_xfr
, sumx
,
1135 nr_ffr
, sumf
, bDim
, scale
, oenv
);
1139 view_all(oenv
, NFILE
, fnm
);
1142 // Free index and isize only if they are distinct from index0 and isize0
1145 for (int i
= 0; i
< ngroups
; i
++)
1152 for (int i
= 0; i
< ngroups
; i
++)
1160 done_filenms(NFILE
, fnm
);
1162 output_env_done(oenv
);