2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Implements gmx::analysismodules::Freevolume.
39 * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
40 * \ingroup module_trajectoryanalysis
49 #include <gromacs/trajectoryanalysis.h>
50 #include <gromacs/pbcutil/pbc.h>
51 #include <gromacs/pbcutil/rmpbc.h>
52 #include <gromacs/utility/smalloc.h>
53 #include <gromacs/math/vectypes.h>
54 #include <gromacs/math/vec.h>
55 #include <gromacs/math/do_fit.h>
57 #define MAX_NTRICVEC 12
69 void make_graph(int mgwi_natoms, rvec *mgwi_x, std::vector< std::vector< node > > &mgwi_graph)
71 mgwi_graph.resize(mgwi_natoms);
72 for (int i = 0; i < mgwi_natoms; i++) {
73 mgwi_graph[i].resize(mgwi_natoms);
75 for (int i = 0; i < mgwi_natoms; i++) {
76 for (int j = 0; j < mgwi_natoms; j++) {
77 rvec_sub(mgwi_x[i], mgwi_x[j], mgwi_graph[i][j].r);
78 mgwi_graph[i][j].n = 0;
83 void update_graph(std::vector< std::vector< node > > &ugwi_graph, rvec *ugwi_x, long double ugwi_epsi) {
85 int ugwi_for = ugwi_graph.size();
86 for (int i = 0; i < ugwi_for; i++) {
87 for (int j = i; j < ugwi_for; j++) {
88 rvec_sub(ugwi_x[i], ugwi_x[j], ugwi_temp);
89 rvec_dec(ugwi_temp, ugwi_graph[i][j].r.as_vec());
90 if (norm(ugwi_temp) <= ugwi_epsi) {
103 void check_domains(long double cd_delta, int cd_frames, std::vector< std::vector< std::vector< node > > > &cd_graph) {
104 int cd_for1 = cd_graph.size(), cd_for2 = cd_graph[1].size();
105 for (int k = 0; k < cd_for1; k++) {
106 for (int i = 0; i < cd_for2; i++) {
107 for (int j = 0; j < cd_for2; j++) {
108 if (cd_graph[k][i][j].n >= cd_frames * cd_delta) {
109 cd_graph[k][i][j].yep = true;
112 cd_graph[k][i][j].yep = false;
119 void find_domain_sizes(std::vector< std::vector< std::vector< node > > > fds_graph, std::vector< std::vector< int > > &fds_domsizes) {
120 fds_domsizes.resize(fds_graph.size());
121 int fds_for1 = fds_graph.size(), fds_for2 = fds_graph[1].size();
122 for (int i = 0; i < fds_for1; i++) {
123 fds_domsizes[i].resize(fds_for2, 0);
124 for (int j = 0; j < fds_for2; j++) {
125 for (int k = 0; k < fds_for2; k++) {
126 if (fds_graph[i][j][k].yep) {
127 fds_domsizes[i][j]++;
134 void get_maxsized_domain(std::vector< int > &gmd_max_d, std::vector< std::vector< std::vector< node > > > gmd_graph, std::vector< std::vector< int > > gmd_domsizes) {
135 int gmd_number1 = 0, gmd_number2 = 0;
136 int gmd_for1 = gmd_domsizes.size(), gmd_for2 = gmd_domsizes[0].size(), gmd_for3 = gmd_graph[1][1].size();
137 for (int i = 0; i < gmd_for1; i++) {
138 for (int j = 0; j < gmd_for2; j++) {
139 if (gmd_domsizes[i][j] >= gmd_domsizes[gmd_number1][gmd_number2]) {
146 for (int i = 0; i < gmd_for3; i++) {
147 if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
148 gmd_max_d.push_back(i);
153 void delete_domain_from_graph(std::vector< std::vector< std::vector< node > > > &ddf_graph, std::vector< int > ddf_domain) {
154 int ddfg_for1 = ddf_domain.size(), ddfg_for2 = ddf_graph.size(), ddfg_for3 = ddf_graph[1].size();
155 for (int i = 0; i < ddfg_for1; i++) {
156 for (int k = 0; k < ddfg_for2; k++) {
157 for (int j = 0; j < ddfg_for3; j++) {
158 if (ddf_graph[k][ddf_domain[i]][j].yep) {
159 ddf_graph[k][ddf_domain[i]][j].yep = false;
161 if (ddf_graph[k][j][ddf_domain[i]].yep) {
162 ddf_graph[k][j][ddf_domain[i]].yep = false;
169 bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain_min_size) {
170 int cd_for1 = cd_domsizes.size(), cd_for2 = cd_domsizes[0].size();
171 for (int i = 0; i < cd_for1; i++) {
172 for (int j = 0; j < cd_for2; j++) {
173 if (cd_domsizes[i][j] >= cd_domain_min_size) {
181 void print_domains(std::vector< std::vector< int > > pd_domains, std::string fnNdx_) {
182 std::freopen(fnNdx_.c_str(), "w+", stdout);
184 for (int i = 0; i < pd_domains.size(); i++) {
185 std::cout << "[domain_" << i + 1 << "]\n";
187 for (int j = 0; j < pd_domains[i].size(); j++) {
189 if (write_count > 20) {
193 std::cout << pd_domains[i][j] << " ";
201 * Class used to compute free volume in a simulations box.
203 * Inherits TrajectoryAnalysisModule and all functions from there.
204 * Does not implement any new functionality.
206 * \ingroup module_trajectoryanalysis
208 class Domains : public TrajectoryAnalysisModule
215 //! Set the options and setting
216 virtual void initOptions(IOptionsContainer *options,
217 TrajectoryAnalysisSettings *settings);
219 //! First routine called by the analysis framework
220 // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
221 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
222 const TopologyInformation &top);
224 virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
225 const t_trxframe &fr);
227 //! Call for each frame of the trajectory
228 // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
229 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
230 TrajectoryAnalysisModuleData *pdata);
232 //! Last routine called by the analysis framework
233 // virtual void finishAnalysis(t_pbc *pbc);
234 virtual void finishAnalysis(int nframes);
236 //! Routine to write output, that is additional over the built-in
237 virtual void writeOutput();
242 std::vector< std::vector< std::vector< node > > > graph;
243 std::vector< std::vector< int > > domains;
244 std::vector< std::vector< int > > domsizes;
245 std::vector< int > index;
246 std::vector< int > numbers;
247 std::vector< std::vector < RVec > > trajectory;
250 int etalon_frame = 0; // should be selectable
251 int domain_min_size = 5; // should be selectable
252 const long double delta = 0.95; //0.95 // should be selectable
253 const long double epsi = 0.1; //0.3 колебания внутри домена // should be selectable
255 // Copy and assign disallowed by base.
258 Domains::Domains(): TrajectoryAnalysisModule()
267 Domains::initOptions(IOptionsContainer *options,
268 TrajectoryAnalysisSettings *settings)
270 static const char *const desc[] = {
271 "[THISMODULE] to be done"
273 // Add the descriptive text (program help text) to the options
274 settings->setHelpText(desc);
275 // Add option for selecting a subset of atoms
276 options->addOption(SelectionOption("select")
277 .store(&selec).required()
278 .description("Atoms that are considered as part of the excluded volume"));
279 // Add option for output file name
280 options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
281 .store(&fnNdx_).defaultBasename("rcore")
282 .description("Index file from the rcore"));
283 // Control input settings
284 settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
285 settings->setPBC(true);
289 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
290 const TopologyInformation &top)
292 domains_ePBC = top.ePBC();
296 Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
297 const t_trxframe &fr)
300 t_pbc *ppbc = settings.hasPBC() ? &pbc : NULL;
302 copy_mat(fr.box, boxx);
304 set_pbc(ppbc, domains_ePBC, boxx);
306 ConstArrayRef< int > atomind = selec.atomIndices();
308 for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
309 index.push_back(*ai);
311 trajectory.resize(frames + 1);
312 trajectory[frames].resize(selec.atomCount());
314 for (int i = 0; i < selec.atomCount(); i++) {
315 trajectory[frames][i] = fr.x[index[i]];
322 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
323 TrajectoryAnalysisModuleData *pdata)
325 trajectory.resize(frames + 1);
326 trajectory[frames].resize(index.size());
327 for (int i = 0; i < index.size(); i++) {
328 trajectory[frames][i] = fr.x[index[i]];
333 //domains -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
334 //domains -s '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.tpr' -f '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.10th.xtc' -select 'name CA'
337 Domains::finishAnalysis(int nframes)
340 int bone = index.size() - domain_min_size + 1;
343 rvec **traj, **etalon;
347 for (int i = 0; i < bone; i++) {
348 snew(w_rls[i], index.size());
349 for (int j = 0; j < index.size(); j++) {
350 if (j >= i && j <= i + domain_min_size - 1) {
359 for (int i = 0; i < bone; i++) {
360 snew(etalon[i], index.size());
361 for (int j = 0; j < index.size(); j++) {
362 copy_rvec(trajectory[etalon_frame][j].as_vec(), etalon[i][j]);
366 for (int i = 0; i < bone; i++) {
367 reset_x(index.size(), NULL, index.size(), NULL, etalon[i], w_rls[i]);
368 make_graph(index.size(), etalon[i], graph[i]);
371 std::chrono::time_point<std::chrono::system_clock> start1, end1;
372 start1 = std::chrono::system_clock::now();
374 for (int i = 1; i < frames; i++) {
375 for (int j = 0; j < bone; j++) {
377 snew(traj[j], index.size());
378 for (int k = 0; k < index.size(); k++) {
379 copy_rvec(trajectory[i][k].as_vec(), traj[j][k]);
384 #pragma omp for schedule(dynamic)
385 for (int j = 0; j < bone; j++) {
386 //#pragma omp ordered
388 reset_x(index.size(), NULL, index.size(), NULL, etalon[j], w_rls[j]);
389 reset_x(index.size(), NULL, index.size(), NULL, traj[j], w_rls[j]);
390 do_fit(index.size(), w_rls[j], etalon[j], traj[j]);
391 update_graph(graph[j], traj[j], epsi);
395 std::cout << i << " / " << frames << " processed\n";
397 end1 = std::chrono::system_clock::now();
398 int seconds = 0, minutes = 0, hours = 0;
399 seconds = std::chrono::duration_cast<std::chrono::seconds>(end1-start1).count();
400 seconds = (float)seconds * (float)( (float)(frames - i) / (float)i );
401 hours = seconds / 3600;
402 seconds -= hours * 3600;
403 minutes = seconds / 60;
405 std::cout << "foretold time: " << hours << " hour(s) " << minutes << " minute(s) " << seconds << "second(s)\n";
408 check_domains(delta, frames, graph);
409 find_domain_sizes(graph, domsizes);
410 std::vector< int > a;
412 while (check_domsizes(domsizes, domain_min_size)) {
413 domains.push_back(a);
414 get_maxsized_domain(domains.back(), graph, domsizes);
415 delete_domain_from_graph(graph, domains.back());
417 find_domain_sizes(graph, domsizes);
419 for (int i = 0; i < bone; i++) {
430 Domains::writeOutput()
432 print_domains(domains, fnNdx_); // see function for details | numbers from index
433 std::cout << "\n END \n";
438 * The main function for the analysis template.
441 main(int argc, char *argv[])
443 return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);