2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/fileio/confio.h"
45 #include "gromacs/fileio/trxio.h"
46 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/gmxana/gmx_ana.h"
48 #include "gromacs/gmxana/gstat.h"
49 #include "gromacs/gmxana/princ.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/pbcutil/rmpbc.h"
52 #include "gromacs/topology/index.h"
53 #include "gromacs/topology/topology.h"
54 #include "gromacs/utility/arraysize.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/futil.h"
57 #include "gromacs/utility/smalloc.h"
60 static void calc_principal_axes(const t_topology
* top
, rvec
* x
, int* index
, int n
, matrix axes
, rvec inertia
)
64 sub_xcm(x
, n
, index
, top
->atoms
.atom
, xcm
, FALSE
);
65 principal_comp(n
, index
, top
->atoms
.atom
, x
, axes
, inertia
);
68 int gmx_principal(int argc
, char* argv
[])
70 const char* desc
[] = {
71 "[THISMODULE] calculates the three principal axes of inertia for a group",
72 "of atoms. NOTE: Old versions of GROMACS wrote the output data in a",
73 "strange transposed way. As of GROMACS 5.0, the output file paxis1.dat",
74 "contains the x/y/z components of the first (major) principal axis for",
75 "each frame, and similarly for the middle and minor axes in paxis2.dat",
78 static gmx_bool foo
= FALSE
;
80 t_pargs pa
[] = { { "-foo", FALSE
, etBOOL
, { &foo
}, "Dummy option to avoid empty array" } };
97 gmx_output_env_t
* oenv
;
98 gmx_rmpbc_t gpbc
= nullptr;
101 t_filenm fnm
[] = { { efTRX
, "-f", nullptr, ffREAD
}, { efTPS
, nullptr, nullptr, ffREAD
},
102 { efNDX
, nullptr, nullptr, ffOPTRD
}, { efXVG
, "-a1", "paxis1", ffWRITE
},
103 { efXVG
, "-a2", "paxis2", ffWRITE
}, { efXVG
, "-a3", "paxis3", ffWRITE
},
104 { efXVG
, "-om", "moi", ffWRITE
} };
105 #define NFILE asize(fnm)
107 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_TIME_UNIT
| PCA_CAN_VIEW
, NFILE
, fnm
,
108 asize(pa
), pa
, asize(desc
), desc
, 0, nullptr, &oenv
))
114 for (i
= 0; i
< DIM
; i
++)
116 snew(legend
[i
], STRLEN
);
117 sprintf(legend
[i
], "%c component", 'X' + i
);
120 axis1
= xvgropen(opt2fn("-a1", NFILE
, fnm
), "Principal axis 1 (major axis)",
121 output_env_get_xvgr_tlabel(oenv
), "Component (nm)", oenv
);
122 xvgr_legend(axis1
, DIM
, legend
, oenv
);
124 axis2
= xvgropen(opt2fn("-a2", NFILE
, fnm
), "Principal axis 2 (middle axis)",
125 output_env_get_xvgr_tlabel(oenv
), "Component (nm)", oenv
);
126 xvgr_legend(axis2
, DIM
, legend
, oenv
);
128 axis3
= xvgropen(opt2fn("-a3", NFILE
, fnm
), "Principal axis 3 (minor axis)",
129 output_env_get_xvgr_tlabel(oenv
), "Component (nm)", oenv
);
130 xvgr_legend(axis3
, DIM
, legend
, oenv
);
132 sprintf(legend
[XX
], "Axis 1 (major)");
133 sprintf(legend
[YY
], "Axis 2 (middle)");
134 sprintf(legend
[ZZ
], "Axis 3 (minor)");
136 fmoi
= xvgropen(opt2fn("-om", NFILE
, fnm
), "Moments of inertia around inertial axes",
137 output_env_get_xvgr_tlabel(oenv
), "I (au nm\\S2\\N)", oenv
);
138 xvgr_legend(fmoi
, DIM
, legend
, oenv
);
140 for (i
= 0; i
< DIM
; i
++)
146 read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), &top
, &ePBC
, nullptr, nullptr, box
, TRUE
);
148 get_index(&top
.atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &gnx
, &index
, &grpname
);
150 natoms
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
152 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, natoms
);
156 gmx_rmpbc(gpbc
, natoms
, box
, x
);
158 calc_principal_axes(&top
, x
, index
, gnx
, axes
, moi
);
160 fprintf(axis1
, "%15.10f %15.10f %15.10f %15.10f\n", t
, axes
[XX
][XX
], axes
[XX
][YY
],
162 fprintf(axis2
, "%15.10f %15.10f %15.10f %15.10f\n", t
, axes
[YY
][XX
], axes
[YY
][YY
],
164 fprintf(axis3
, "%15.10f %15.10f %15.10f %15.10f\n", t
, axes
[ZZ
][XX
], axes
[ZZ
][YY
],
166 fprintf(fmoi
, "%15.10f %15.10f %15.10f %15.10f\n", t
, moi
[XX
], moi
[YY
], moi
[ZZ
]);
167 } while (read_next_x(oenv
, status
, &t
, x
, box
));
169 gmx_rmpbc_done(gpbc
);