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
51 #include <gromacs/trajectoryanalysis.h>
52 #include <gromacs/pbcutil/pbc.h>
53 #include <gromacs/utility/smalloc.h>
54 #include <gromacs/math/do_fit.h>
55 #include <gromacs/fileio/trxio.h>
61 double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
62 return sqrt( pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
63 pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
64 pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2) );
67 void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
68 double t01, t02, t03, t04, t05, t06, t07;
69 t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
70 t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
71 t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
72 t04 = sqrt(t01 + t02 + t03);
73 t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
74 t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
75 t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
77 //#pragma omp parallel sections
82 Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
84 Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
86 Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
88 Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
89 2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
90 2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
92 Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
93 2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
94 2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
96 Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
97 2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
101 /*void searchL (double &fLd, double &fLdd, double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C, double L, double fx, double fy, double fz, double fa, double fb, double fc) {
102 double AmLfa, BmLfb, CmLfc, BIXpXmLfx, BIYpYmLfy, BIZpZmLfz, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20;
106 BIXpXmLfx = bix + x - L * fx;
107 BIYpYmLfy = biy + y - L * fy;
108 BIZpZmLfz = biz + z - L * fz;
109 t1 = (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc));
110 t2 = (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc));
111 t3 = (cos(AmLfa) * sin(CmLfc) - cos(CmLfc) * sin(AmLfa) * sin(BmLfb));
112 t4 = (sin(AmLfa) * sin(CmLfc) + cos(AmLfa) * cos(CmLfc) * sin(BmLfb));
113 t5 = pow(aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy, 2);
114 t6 = pow(aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx, 2);
115 t7 = pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2);
117 t8 = fc * sin(AmLfa) * sin(CmLfc);
118 t9 = fa * cos(AmLfa) * cos(CmLfc);
119 t10 = fb * cos(AmLfa) * cos(BmLfb) * sin(CmLfc);
120 t11 = fc * cos(AmLfa) * cos(CmLfc) * sin(BmLfb);
121 t12 = fa * sin(AmLfa) * sin(BmLfb) * sin(CmLfc);
123 t13 = fa * cos(AmLfa) * sin(BmLfb) * sin(CmLfc);
124 t14 = fc * cos(AmLfa) * sin(CmLfc);
125 t15 = fa * cos(CmLfc) * sin(AmLfa);
126 t16 = fb * cos(BmLfb) * sin(AmLfa) * sin(CmLfc);
127 t17 = fc * cos(CmLfc) * sin(AmLfa) * sin(BmLfb);
129 t18 = fx * cos(BmLfb) * sin(CmLfc);
130 t19 = fc * cos(BmLfb) * cos(CmLfc);
131 t20 = fb * sin(BmLfb) * sin(CmLfc);
133 fLd += (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
134 ( BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
135 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
136 (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) -
137 fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy + fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
138 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
139 (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
140 BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
141 fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * t3 +
142 fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) /
143 (2 * sqrt( t5 + t6 + pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2)));
146 fLdd += ( (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
147 (2 * fb * fx * cos(BmLfb) - pow(fb, 2) * sin(BmLfb) * BIXpXmLfx - 2 * fa * fy * cos(AmLfa) * cos(BmLfb) + 2 * fa * fz * cos(BmLfb) * sin(AmLfa) + 2 * fb * fz * cos(AmLfa) * sin(BmLfb) +
148 2 * fb * fy * sin(AmLfa) * sin(BmLfb) + pow(fa, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz + pow(fb, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz +
149 pow(fa, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + pow(fb, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + 2 * fa * fb * cos(AmLfa) * sin(BmLfb) * BIYpYmLfy -
150 2 * fa * fb * sin(AmLfa) * sin(BmLfb) * BIZpZmLfz) +
151 (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
152 fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
153 BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
154 fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * (cos(AmLfa) * sin(CmLfc) -
155 cos(CmLfc) * sin(AmLfa) * sin(BmLfb)) + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) -
156 fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
157 (2 * BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) +
158 fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) -
159 fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) +
160 t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
161 2 * fy * t3 + 2 * fz * t4 +
162 2 * fx * cos(BmLfb) * cos(CmLfc) - 2 * fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - 2 * fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
163 (2 * aix + 2 * t3 * BIYpYmLfy - 2 * (sin(AmLfa) * sin(CmLfc) +
164 cos(AmLfa) * cos(CmLfc) * sin(BmLfb)) * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
165 (BIZpZmLfz * (pow(fa, 2) * sin(AmLfa) * sin(CmLfc) +
166 pow(fc, 2) * sin(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(AmLfa) * cos(CmLfc) + pow(fa, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) +
167 pow(fb, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + 2 * fa * fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) +
168 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * sin(CmLfc) - 2 * fa * fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + BIYpYmLfy * (pow(fa, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) -
169 pow(fc, 2) * cos(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(CmLfc) * sin(AmLfa) - pow(fa, 2) * cos(AmLfa) * sin(CmLfc) + pow(fb, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) +
170 pow(fc, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) - 2 * fa * fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
171 2 * fb * fc * cos(BmLfb) * sin(AmLfa) * sin(CmLfc)) - 2 * fz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
172 fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
173 2 * fy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) -
174 fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * fb * fx * cos(CmLfc) * sin(BmLfb) + 2 * fc * t18 + pow(fb, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx +
175 pow(fc, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx - 2 * fb * fc * sin(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
176 (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) *
177 (2 * BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + 2 * fy * t1 - 2 * fz * t2 + 2 * BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + 2 * t18 + 2 * t19 * BIXpXmLfx - 2 * t20 * BIXpXmLfx) +
178 (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
179 fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) *
180 (2 * fx * sin(BmLfb) - 2 * fy * cos(BmLfb) * sin(AmLfa) + 2 * fb * cos(BmLfb) * BIXpXmLfx - 2 * fz * cos(AmLfa) * cos(BmLfb) - 2 * fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
181 2 * fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + 2 * fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + 2 * fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
182 (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
183 (BIZpZmLfz * (pow(fa, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) -
184 pow(fc, 2) * cos(CmLfc) * sin(AmLfa) - 2 * fa * t14 - pow(fa, 2) * cos(CmLfc) * sin(AmLfa) + pow(fb, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
185 pow(fc, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) - 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * t16 +
186 2 * fa * t17) + BIYpYmLfy * (pow(fa, 2) * cos(AmLfa) * cos(CmLfc) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) - 2 * fa * t8 +
187 pow(fa, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fb, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fc, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) -
188 2 * fa * t10 - 2 * fa * t11 - 2 * fb * t19 * sin(AmLfa)) - 2 * fz * (t8 - t9 + t10 + t11 - t12) - 2 * fy * (t13 - t14 - t15 + t16 + t17) -
189 2 * fc * fx * cos(BmLfb) * cos(CmLfc) + 2 * fb * fx * sin(BmLfb) * sin(CmLfc) +
190 pow(fb, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + pow(fc, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + 2 * fb * fc * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx)) /
191 (2 * sqrt( t5 + t6 + t7)) -
192 ( (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
193 (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
194 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
195 (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
196 fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
197 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
198 (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
199 BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
200 fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) *
201 ( (2 * aix + 2 * t3 * BIYpYmLfy - 2 * t4 * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
202 (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
203 BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
204 fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
205 (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * t2 * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
206 (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) +
207 t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) - (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
208 (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
209 fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy))) /
210 (4 * pow ( t5 + t6 + t7, 3 / 2));
214 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
217 #pragma omp for schedule(dynamic)
218 for (int i = 0; i < b.size(); i++) {
220 b[i][0] = (temp[2] + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (temp[1] + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (temp[0] + x);
221 b[i][1] = (temp[1] + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (temp[2] + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (temp[0] + x);
222 b[i][2] = cos(A) * cos(B) * (temp[2] + z) - sin(B) * (temp[0] + x) + cos(B) * sin(A) * (temp[1] + y);
227 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
236 for (int i = 0; i < pairs.size(); i++) {
237 rvec_inc(midA, a[pairs[i].first]);
238 rvec_inc(midB, b[pairs[i].second]);
240 midA[0] /= pairs.size();
241 midA[1] /= pairs.size();
242 midA[2] /= pairs.size();
244 midB[0] /= pairs.size();
245 midB[1] /= pairs.size();
246 midB[2] /= pairs.size();
249 /*float*/double FrameMidLength(std::vector< RVec > a) {
254 for (int i = 0; i < a.size(); i++) {
260 void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs, double FtCnst) {
261 double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
263 CalcMid(a, b, ma, mb, pairs);
265 double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
266 for (int i = 0; i < pairs.size(); i++) {
267 f1 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
276 while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
277 f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
278 for (int i = 0; i < pairs.size(); i++) {
279 searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
283 for (int i = 0; i < pairs.size(); i++) {
284 f2 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + (x - l * fx), b[pairs[i].second][1] + (y - l * fy), b[pairs[i].second][2] + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
288 x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
289 fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
292 l *= 0.99; // may be put somewhere l > 0
296 //std::cout << " " << count << "\n";
297 ApplyFit(b, x, y, z, A, B, C);
301 void printPDBtraj(const char* Fname, std::vector< std::vector < RVec > > trj, std::vector< std::vector< std::pair< int, int > > > prs, std::vector< int > ndx) {
302 FILE* a = fopen (Fname, "a");
303 for (int i = 1; i < trj.size(); i++) {
304 fprintf(a, "TITLE test t= %9.5f step= %5d\nMODEL %d\n", (/*float*/double)(i + 1), i + 1, i + 1);
305 for (int j = 0; j < prs.size(); j++) {
306 for (int k = 0; k < prs[j].size(); k++) {
307 fprintf(a, "ATOM %5d CA CA %5d %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s \n", ndx[prs[j][k].first], ndx[prs[j][k].first], (trj[i][prs[j][k].first][0] * 10), (trj[i][prs[j][k].first][1] * 10), (trj[i][prs[j][k].first][2] * 10), 1.0, 20.0, " ", " ");
310 fprintf(a, "TER\nENDMDL\n");
312 //fprintf(a, "END\n");
316 bool myfunction (const std::vector< std::pair< int, int > > i, const std::vector< std::pair< int, int > > j) {
317 return i.size() < j.size();
320 void domain_reading(SelectionList sList, std::vector< int > ind, std::vector< std::vector < std::pair< int, int > > > &d_ind_num) {
321 d_ind_num.resize(sList.size() - 1);
322 for (int i = 1; i < sList.size(); i++) {
323 d_ind_num[i - 1].resize(0);
324 for (ArrayRef< const int >::iterator ai = sList[i].atomIndices().begin(); (ai < sList[i].atomIndices().end()); ai++) {
325 d_ind_num[i - 1].push_back(std::make_pair(*ai, 0));
329 for (int i = 0; i < d_ind_num.size(); i++) {
331 for (int j = 0; j < d_ind_num[i].size(); j++) {
333 while (ind[k] != d_ind_num[i][j].first) {
336 d_ind_num[i][j].second = k;
340 // its a point of interest as we can either take biggest domains of the given set and work from those or we take them in given order
341 //std::sort(d_ind_num.begin(), d_ind_num.end(), myfunction);
343 for (int i = 0; i < d_ind_num.size(); i++) {
344 if (d_ind_num[i].size() >= 4) {
345 for (int k = 0; k < d_ind_num[i].size(); k++) {
346 for (int j = i + 1; j < d_ind_num.size(); j++) {
347 for (int x = 0; x < d_ind_num[j].size(); x++) {
348 if (d_ind_num[j][x].first == d_ind_num[i][k].first) {
349 d_ind_num[j].erase(d_ind_num[j].begin() + x);
359 void domain_expansion(std::vector< int > indx, std::vector< std::vector < std::pair< int, int > > > &d_ind_num, std::vector< RVec > frame) {
360 std::vector< bool > flag;
361 flag.resize(indx.size(), 1);
362 for (int i = 0; i < d_ind_num.size(); i++) {
363 for (int j = 0; j < d_ind_num[i].size(); j++) {
364 flag[d_ind_num[i][j].second] = 0;
367 /*float*/double dist;
369 for (int i = 0; i < indx.size(); i++) {
372 for (int x1 = 0; x1 < d_ind_num.size(); x1++) {
373 for (int x2 = 0; x2 < d_ind_num[x1].size(); x2++) {
375 dist = (frame[i] - frame[d_ind_num[0][0].first]).norm();
378 if ((frame[i] - frame[d_ind_num[x1][x2].first]).norm() < dist) {
380 dist = (frame[i] - frame[d_ind_num[x1][x2].first]).norm();
384 d_ind_num[x3].push_back(std::make_pair(indx[i], i));
389 void MakeDomainFitPairs(std::vector< std::vector< std::pair< int, int > > > &p, std::vector< std::vector< std::pair< int, int > > > d) {
391 for (int i = 0; i < d.size(); i++) {
393 for (int j = 0; j < d[i].size(); j++) {
394 p[i].push_back(std::make_pair(d[i][j].second, d[i][j].second));
399 void TrajFitting(std::vector< std::vector < RVec > > &trj, double FC, std::vector< std::vector< std::pair< int, int > > > prs) {
401 // its needed to be thought through on whats shared and whats private
402 //#pragma omp parallel for shared(trj, prs, FC)/* private(clone)*/
403 for (int i = 1; i < trj.size()/*2*/; i++) {
404 std::vector< std::vector < RVec > > clone;
406 clone.resize(prs.size(), trj[i]);
407 //std::chrono::time_point<std::chrono::system_clock> now1, now2;
408 for (int j = 0; j < prs.size(); j++) {
409 //now1 = std::chrono::system_clock::now();
410 MyFitNew(trj[0], clone[j], prs[j], FC);
411 //now2 = std::chrono::system_clock::now();
412 //std::cout << " " << std::chrono::duration_cast<std::chrono::milliseconds>(now2 - now1).count() << "\n";
413 for (int k = 0; k < prs[j].size(); k++) {
414 trj[i][prs[j][k].first] = clone[j][prs[j][k].first];
418 //std::cout << "frame " << i << " fitted\n";
422 class Fitng : public TrajectoryAnalysisModule
429 //! Set the options and setting
430 virtual void initOptions(IOptionsContainer *options,
431 TrajectoryAnalysisSettings *settings);
433 //! First routine called by the analysis framework
434 // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
435 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
436 const TopologyInformation &top);
438 virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
439 const t_trxframe &fr);
441 //! Call for each frame of the trajectory
442 // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
443 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
444 TrajectoryAnalysisModuleData *pdata);
446 //! Last routine called by the analysis framework
447 // virtual void finishAnalysis(t_pbc *pbc);
448 virtual void finishAnalysis(int nframes);
450 //! Routine to write output, that is additional over the built-in
451 virtual void writeOutput();
456 std::vector< int > index;
458 std::vector< std::vector< std::pair< int, int > > > domains_index_number;
459 std::vector< std::vector < RVec > > trajectory;
460 std::vector< t_trxframe > saved_traj;
462 std::vector< std::vector< std::pair< int, int > > > pairs;
466 std::string OutPutTrjName;
470 Fitng::Fitng(): TrajectoryAnalysisModule()
479 Fitng::initOptions(IOptionsContainer *options,
480 TrajectoryAnalysisSettings *settings)
482 static const char *const desc[] = {
483 "[THISMODULE] to be done"
485 // Add the descriptive text (program help text) to the options
486 settings->setHelpText(desc);
487 // Add option for selecting a subset of atoms
488 //options->addOption(SelectionOption("select")
489 // .store(&selec).required()
490 // .description("Atoms that are considered as part of the excluded volume"));
491 // Add option for selection list
492 options->addOption(SelectionOption("select_residue_and_domains").storeVector(&sel_)
493 .required().dynamicMask().multiValue()
494 .description("Selected group and domains based on it"));
495 // Add option for Fit constant
496 options->addOption(DoubleOption("FitConst")
498 .description("Fitting untill diff <= FitConst, by default == 0.00001"));
499 // Add option for output pdb traj (meh edition)
500 options->addOption(StringOption("pdb_out")
501 .store(&OutPutTrjName)
502 .description("transformed trajectory"));
503 // Add option for trajectory Fit transformation
504 options->addOption(IntegerOption("method")
506 .description("1 - fitting all -> all | 2 - keep only domains, all -> all | 3 - keep only domains, domain -> domain | 4 - extend domains to full coverage, domain -> domain"));
507 // Control input settings
508 settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
509 settings->setPBC(true);
512 // /home/toluk/Develop/samples/reca_test_ground/
513 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL0' -n '/home/toluk/Develop/samples/reca_test_ground/t0.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/FullCAFit.pdb' -e 1000
514 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsFit.pdb' -e 1000 -method 2
515 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsOnlyFit.pdb' -e 1000 -method 3
516 // -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsWExtraOnlyFit.pdb' -e 1000 -method 4
518 // -s '/home/titov_ai/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/titov_ai/Develop/samples/reca_rd/reca_rd.mono.xtc' -sf '/home/titov_ai/Develop/samples/reca_rd/SLCa' -n '/home/titov_ai/Develop/samples/reca_rd/CA.ndx' -pdb_out '/home/titov_ai/Develop/samples/reca_rd/test.xtc' -method 1
520 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/test.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/testfit.xtc' -method 1
523 Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
526 for (ArrayRef< const int >::iterator ai = sel_[0].atomIndices().begin(); (ai < sel_[0].atomIndices().end()); ai++) {
527 index.push_back(*ai);
530 saved_traj.resize(0);
532 if (sel_.size() > 1 && method >= 2) {
533 domain_reading(sel_, index, domains_index_number);
539 Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
541 trajectory.resize(0);
542 std::vector < RVec > a;
544 trajectory.push_back(a);
545 for (int i = 0; i < index.size(); i++) {
546 trajectory[0].push_back(fr.x[index[i]]);
548 trajectory.push_back(a);
552 for (int i = 0; i < index.size(); i++) {
553 pairs[0].push_back(std::make_pair(i, i));
558 for (int i = 0; i < domains_index_number.size(); i++) {
559 for (int j = 0; j < domains_index_number[i].size(); j++) {
560 pairs[0].push_back(std::make_pair(domains_index_number[i][j].second, domains_index_number[i][j].second));
565 MakeDomainFitPairs(pairs, domains_index_number);
568 domain_expansion(index, domains_index_number, trajectory[0]);
569 MakeDomainFitPairs(pairs, domains_index_number);
573 op = open_trx(OutPutTrjName.c_str(), "w+");
577 Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
578 TrajectoryAnalysisModuleData *pdata)
580 trajectory[1].resize(0);
581 for (int i = 0; i < index.size(); i++) {
582 trajectory[1].push_back(fr.x[index[i]]);
584 TrajFitting(trajectory, FitConst, pairs);
587 for (int i = 0; i < index.size(); i++) {
588 copy_rvec(trajectory[1][i].as_vec(), c.x[index[i]]);
590 //trajectory.push_back(trajectory[1]);
591 if (write_trxframe_indexed(op, static_cast<const t_trxframe *>(&c), index.size(), static_cast<const int *>(&index.front())/*saved_traj[i].index*/, NULL) != 0) {
592 std::cout << "\nFileOutPut error\nframe number = " << frnr << "\n";
595 printPDBtraj((OutPutTrjName + ".pdb").c_str(), trajectory, pairs, index);
600 Fitng::finishAnalysis(int nframes)
602 std::cout << "\nGame Over\n";
608 //printPDBtraj((OutPutTrjName + ".pdb").c_str(), trajectory, pairs, index);
613 * The main function for the analysis template.
616 main(int argc, char *argv[])
618 return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Fitng>(argc, argv);