?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Cleaned up memory usage in gmx traj and trjconv
[gromacs.git] / src / gromacs / gmxana / gmx_helixorient.cpp
blob8b5a3c9bcfd1135bb7b260326abd3f86f15022cc
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,2015,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.
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 #include "gmxpre.h"
38
39 #include <cmath>
40
41 #include "gromacs/commandline/pargs.h"
42 #include "gromacs/fileio/trxio.h"
43 #include "gromacs/fileio/xvgr.h"
44 #include "gromacs/gmxana/gmx_ana.h"
45 #include "gromacs/gmxana/gstat.h"
46 #include "gromacs/math/do_fit.h"
47 #include "gromacs/math/utilities.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/pbcutil/pbc.h"
50 #include "gromacs/pbcutil/rmpbc.h"
51 #include "gromacs/topology/index.h"
52 #include "gromacs/topology/topology.h"
53 #include "gromacs/utility/arraysize.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/smalloc.h"
57
58 int gmx_helixorient(int argc, char *argv[])
59 {
60 const char *desc[] = {
61 "[THISMODULE] calculates the coordinates and direction of the average",
62 "axis inside an alpha helix, and the direction/vectors of both the",
63 "C[GRK]alpha[grk] and (optionally) a sidechain atom relative to the axis.[PAR]",
64 "As input, you need to specify an index group with C[GRK]alpha[grk] atoms",
65 "corresponding to an [GRK]alpha[grk]-helix of continuous residues. Sidechain",
66 "directions require a second index group of the same size, containing",
67 "the heavy atom in each residue that should represent the sidechain.[PAR]",
68 "[BB]Note[bb] that this program does not do any fitting of structures.[PAR]",
69 "We need four C[GRK]alpha[grk] coordinates to define the local direction of the helix",
70 "axis.[PAR]",
71 "The tilt/rotation is calculated from Euler rotations, where we define",
72 "the helix axis as the local [IT]x[it]-axis, the residues/C[GRK]alpha[grk] vector as [IT]y[it], and the",
73 "[IT]z[it]-axis from their cross product. We use the Euler Y-Z-X rotation, meaning",
74 "we first tilt the helix axis (1) around and (2) orthogonal to the residues",
75 "vector, and finally apply the (3) rotation around it. For debugging or other",
76 "purposes, we also write out the actual Euler rotation angles as [TT]theta[1-3].xvg[tt]"
77 };
78
79 t_topology *top = nullptr;
80 real t;
81 rvec *x = nullptr;
82 matrix box;
83 t_trxstatus *status;
84 int natoms;
85 real theta1, theta2, theta3;
86
87 int i, j, teller = 0;
88 int iCA, iSC;
89 int *ind_CA;
90 int *ind_SC;
91 char *gn_CA;
92 char *gn_SC;
93 rvec v1, v2;
94 rvec *x_CA, *x_SC;
95 rvec *r12;
96 rvec *r23;
97 rvec *r34;
98 rvec *diff13;
99 rvec *diff24;
100 rvec *helixaxis;
101 rvec *residuehelixaxis;
102 rvec *residueorigin;
103 rvec *residuevector;
104 rvec *sidechainvector;
105
106 rvec *residuehelixaxis_t0;
107 rvec *residuevector_t0;
108 rvec *axis3_t0;
109 rvec *residuehelixaxis_tlast;
110 rvec *residuevector_tlast;
111 rvec *axis3_tlast;
112 rvec refaxes[3], newaxes[3];
113 rvec unitaxes[3];
114 rvec rot_refaxes[3], rot_newaxes[3];
115
116 real tilt, rotation;
117 rvec *axis3;
118 real *twist, *residuetwist;
119 real *radius, *residueradius;
120 real *rise, *residuerise;
121 real *residuebending;
122
123 real tmp;
124 real weight[3];
125 t_pbc pbc;
126 matrix A;
127
128 FILE *fpaxis, *fpcenter, *fptilt, *fprotation;
129 FILE *fpradius, *fprise, *fptwist;
130 FILE *fptheta1, *fptheta2, *fptheta3;
131 FILE *fpbending;
132 int ePBC;
133
134 gmx_output_env_t *oenv;
135 gmx_rmpbc_t gpbc = nullptr;
136
137 static gmx_bool bSC = FALSE;
138 static gmx_bool bIncremental = FALSE;
139
140 static t_pargs pa[] = {
141 { "-sidechain", FALSE, etBOOL, {&bSC},
142 "Calculate sidechain directions relative to helix axis too." },
143 { "-incremental", FALSE, etBOOL, {&bIncremental},
144 "Calculate incremental rather than total rotation/tilt." },
145 };
146 #define NPA asize(pa)
147
148 t_filenm fnm[] = {
149 { efTPR, nullptr, nullptr, ffREAD },
150 { efTRX, "-f", nullptr, ffREAD },
151 { efNDX, nullptr, nullptr, ffOPTRD },
152 { efDAT, "-oaxis", "helixaxis", ffWRITE },
153 { efDAT, "-ocenter", "center", ffWRITE },
154 { efXVG, "-orise", "rise", ffWRITE },
155 { efXVG, "-oradius", "radius", ffWRITE },
156 { efXVG, "-otwist", "twist", ffWRITE },
157 { efXVG, "-obending", "bending", ffWRITE },
158 { efXVG, "-otilt", "tilt", ffWRITE },
159 { efXVG, "-orot", "rotation", ffWRITE }
160 };
161 #define NFILE asize(fnm)
162
163 if (!parse_common_args(&argc, argv, PCA_CAN_TIME,
164 NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
165 {
166 return 0;
167 }
168
169 top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);
170
171 for (i = 0; i < 3; i++)
172 {
173 weight[i] = 1.0;
174 }
175
176 /* read index files */
177 printf("Select a group of Calpha atoms corresponding to a single continuous helix:\n");
178 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iCA, &ind_CA, &gn_CA);
179 snew(x_CA, iCA);
180 snew(x_SC, iCA); /* sic! */
181
182 snew(r12, iCA-3);
183 snew(r23, iCA-3);
184 snew(r34, iCA-3);
185 snew(diff13, iCA-3);
186 snew(diff24, iCA-3);
187 snew(helixaxis, iCA-3);
188 snew(twist, iCA);
189 snew(residuetwist, iCA);
190 snew(radius, iCA);
191 snew(residueradius, iCA);
192 snew(rise, iCA);
193 snew(residuerise, iCA);
194 snew(residueorigin, iCA);
195 snew(residuehelixaxis, iCA);
196 snew(residuevector, iCA);
197 snew(sidechainvector, iCA);
198 snew(residuebending, iCA);
199 snew(residuehelixaxis_t0, iCA);
200 snew(residuevector_t0, iCA);
201 snew(axis3_t0, iCA);
202 snew(residuehelixaxis_tlast, iCA);
203 snew(residuevector_tlast, iCA);
204 snew(axis3_tlast, iCA);
205 snew(axis3, iCA);
206
207 if (bSC)
208 {
209 printf("Select a group of atoms defining the sidechain direction (1/residue):\n");
210 get_index(&(top->atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &iSC, &ind_SC, &gn_SC);
211 if (iSC != iCA)
212 {
213 gmx_fatal(FARGS, "Number of sidechain atoms (%d) != number of CA atoms (%d)", iSC, iCA);
214 }
215
216 }
217
218 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
219
220 fpaxis = gmx_ffopen(opt2fn("-oaxis", NFILE, fnm), "w");
221 fpcenter = gmx_ffopen(opt2fn("-ocenter", NFILE, fnm), "w");
222 fprise = gmx_ffopen(opt2fn("-orise", NFILE, fnm), "w");
223 fpradius = gmx_ffopen(opt2fn("-oradius", NFILE, fnm), "w");
224 fptwist = gmx_ffopen(opt2fn("-otwist", NFILE, fnm), "w");
225 fpbending = gmx_ffopen(opt2fn("-obending", NFILE, fnm), "w");
226
227 fptheta1 = gmx_ffopen("theta1.xvg", "w");
228 fptheta2 = gmx_ffopen("theta2.xvg", "w");
229 fptheta3 = gmx_ffopen("theta3.xvg", "w");
230
231 if (bIncremental)
232 {
233 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
234 "Incremental local helix tilt", "Time(ps)", "Tilt (degrees)",
235 oenv);
236 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
237 "Incremental local helix rotation", "Time(ps)",
238 "Rotation (degrees)", oenv);
239 }
240 else
241 {
242 fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
243 "Cumulative local helix tilt", "Time(ps)", "Tilt (degrees)", oenv);
244 fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
245 "Cumulative local helix rotation", "Time(ps)",
246 "Rotation (degrees)", oenv);
247 }
248
249 clear_rvecs(3, unitaxes);
250 unitaxes[0][0] = 1;
251 unitaxes[1][1] = 1;
252 unitaxes[2][2] = 1;
253
254 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
255
256 do
257 {
258 /* initialisation for correct distance calculations */
259 set_pbc(&pbc, ePBC, box);
260 /* make molecules whole again */
261 gmx_rmpbc(gpbc, natoms, box, x);
262
263 /* copy coords to our smaller arrays */
264 for (i = 0; i < iCA; i++)
265 {
266 copy_rvec(x[ind_CA[i]], x_CA[i]);
267 if (bSC)
268 {
269 copy_rvec(x[ind_SC[i]], x_SC[i]);
270 }
271 }
272
273 for (i = 0; i < iCA-3; i++)
274 {
275 rvec_sub(x_CA[i+1], x_CA[i], r12[i]);
276 rvec_sub(x_CA[i+2], x_CA[i+1], r23[i]);
277 rvec_sub(x_CA[i+3], x_CA[i+2], r34[i]);
278 rvec_sub(r12[i], r23[i], diff13[i]);
279 rvec_sub(r23[i], r34[i], diff24[i]);
280 /* calculate helix axis */
281 cprod(diff13[i], diff24[i], helixaxis[i]);
282 svmul(1.0/norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
283
284 tmp = cos_angle(diff13[i], diff24[i]);
285 twist[i] = 180.0/M_PI * std::acos( tmp );
286 radius[i] = std::sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
287 rise[i] = std::abs(iprod(r23[i], helixaxis[i]));
288
289 svmul(radius[i]/norm(diff13[i]), diff13[i], v1);
290 svmul(radius[i]/norm(diff24[i]), diff24[i], v2);
291
292 rvec_sub(x_CA[i+1], v1, residueorigin[i+1]);
293 rvec_sub(x_CA[i+2], v2, residueorigin[i+2]);
294 }
295 residueradius[0] = residuetwist[0] = residuerise[0] = 0;
296
297 residueradius[1] = radius[0];
298 residuetwist[1] = twist[0];
299 residuerise[1] = rise[0];
300
301 residuebending[0] = residuebending[1] = 0;
302 for (i = 2; i < iCA-2; i++)
303 {
304 residueradius[i] = 0.5*(radius[i-2]+radius[i-1]);
305 residuetwist[i] = 0.5*(twist[i-2]+twist[i-1]);
306 residuerise[i] = 0.5*(rise[i-2]+rise[i-1]);
307 residuebending[i] = 180.0/M_PI*std::acos( cos_angle(helixaxis[i-2], helixaxis[i-1]) );
308 }
309 residueradius[iCA-2] = radius[iCA-4];
310 residuetwist[iCA-2] = twist[iCA-4];
311 residuerise[iCA-2] = rise[iCA-4];
312 residueradius[iCA-1] = residuetwist[iCA-1] = residuerise[iCA-1] = 0;
313 residuebending[iCA-2] = residuebending[iCA-1] = 0;
314
315 clear_rvec(residueorigin[0]);
316 clear_rvec(residueorigin[iCA-1]);
317
318 /* average helix axes to define them on the residues.
319 * Just extrapolate second first/list atom.
320 */
321 copy_rvec(helixaxis[0], residuehelixaxis[0]);
322 copy_rvec(helixaxis[0], residuehelixaxis[1]);
323
324 for (i = 2; i < iCA-2; i++)
325 {
326 rvec_add(helixaxis[i-2], helixaxis[i-1], residuehelixaxis[i]);
327 svmul(0.5, residuehelixaxis[i], residuehelixaxis[i]);
328 }
329 copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-2]);
330 copy_rvec(helixaxis[iCA-4], residuehelixaxis[iCA-1]);
331
332 /* Normalize the axis */
333 for (i = 0; i < iCA; i++)
334 {
335 svmul(1.0/norm(residuehelixaxis[i]), residuehelixaxis[i], residuehelixaxis[i]);
336 }
337
338 /* calculate vector from origin to residue CA */
339 fprintf(fpaxis, "%15.12g ", t);
340 fprintf(fpcenter, "%15.12g ", t);
341 fprintf(fprise, "%15.12g ", t);
342 fprintf(fpradius, "%15.12g ", t);
343 fprintf(fptwist, "%15.12g ", t);
344 fprintf(fpbending, "%15.12g ", t);
345
346 for (i = 0; i < iCA; i++)
347 {
348 if (i == 0 || i == iCA-1)
349 {
350 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
351 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", 0.0, 0.0, 0.0);
352 fprintf(fprise, "%15.12g ", 0.0);
353 fprintf(fpradius, "%15.12g ", 0.0);
354 fprintf(fptwist, "%15.12g ", 0.0);
355 fprintf(fpbending, "%15.12g ", 0.0);
356 }
357 else
358 {
359 rvec_sub( bSC ? x_SC[i] : x_CA[i], residueorigin[i], residuevector[i]);
360 svmul(1.0/norm(residuevector[i]), residuevector[i], residuevector[i]);
361 cprod(residuehelixaxis[i], residuevector[i], axis3[i]);
362 fprintf(fpaxis, "%15.12g %15.12g %15.12g ", residuehelixaxis[i][0], residuehelixaxis[i][1], residuehelixaxis[i][2]);
363 fprintf(fpcenter, "%15.12g %15.12g %15.12g ", residueorigin[i][0], residueorigin[i][1], residueorigin[i][2]);
364
365 fprintf(fprise, "%15.12g ", residuerise[i]);
366 fprintf(fpradius, "%15.12g ", residueradius[i]);
367 fprintf(fptwist, "%15.12g ", residuetwist[i]);
368 fprintf(fpbending, "%15.12g ", residuebending[i]);
369 }
370 }
371 fprintf(fprise, "\n");
372 fprintf(fpradius, "\n");
373 fprintf(fpaxis, "\n");
374 fprintf(fpcenter, "\n");
375 fprintf(fptwist, "\n");
376 fprintf(fpbending, "\n");
377
378 if (teller == 0)
379 {
380 for (i = 0; i < iCA; i++)
381 {
382 copy_rvec(residuehelixaxis[i], residuehelixaxis_t0[i]);
383 copy_rvec(residuevector[i], residuevector_t0[i]);
384 copy_rvec(axis3[i], axis3_t0[i]);
385 }
386 }
387 else
388 {
389 fprintf(fptilt, "%15.12g ", t);
390 fprintf(fprotation, "%15.12g ", t);
391 fprintf(fptheta1, "%15.12g ", t);
392 fprintf(fptheta2, "%15.12g ", t);
393 fprintf(fptheta3, "%15.12g ", t);
394
395 for (i = 0; i < iCA; i++)
396 {
397 if (i == 0 || i == iCA-1)
398 {
399 tilt = rotation = 0;
400 }
401 else
402 {
403 if (!bIncremental)
404 {
405 /* Total rotation & tilt */
406 copy_rvec(residuehelixaxis_t0[i], refaxes[0]);
407 copy_rvec(residuevector_t0[i], refaxes[1]);
408 copy_rvec(axis3_t0[i], refaxes[2]);
409 }
410 else
411 {
412 /* Rotation/tilt since last step */
413 copy_rvec(residuehelixaxis_tlast[i], refaxes[0]);
414 copy_rvec(residuevector_tlast[i], refaxes[1]);
415 copy_rvec(axis3_tlast[i], refaxes[2]);
416 }
417 copy_rvec(residuehelixaxis[i], newaxes[0]);
418 copy_rvec(residuevector[i], newaxes[1]);
419 copy_rvec(axis3[i], newaxes[2]);
420
421 /* rotate reference frame onto unit axes */
422 calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
423 for (j = 0; j < 3; j++)
424 {
425 mvmul(A, refaxes[j], rot_refaxes[j]);
426 mvmul(A, newaxes[j], rot_newaxes[j]);
427 }
428
429 /* Determine local rotation matrix A */
430 calc_fit_R(3, 3, weight, rot_newaxes, rot_refaxes, A);
431 /* Calculate euler angles, from rotation order y-z-x, where
432 * x is helixaxis, y residuevector, and z axis3.
433 *
434 * A contains rotation column vectors.
435 */
436
437 theta1 = 180.0/M_PI*std::atan2(A[0][2], A[0][0]);
438 theta2 = 180.0/M_PI*std::asin(-A[0][1]);
439 theta3 = 180.0/M_PI*std::atan2(A[2][1], A[1][1]);
440
441 tilt = std::sqrt(theta1*theta1+theta2*theta2);
442 rotation = theta3;
443 fprintf(fptheta1, "%15.12g ", theta1);
444 fprintf(fptheta2, "%15.12g ", theta2);
445 fprintf(fptheta3, "%15.12g ", theta3);
446
447 }
448 fprintf(fptilt, "%15.12g ", tilt);
449 fprintf(fprotation, "%15.12g ", rotation);
450 }
451 fprintf(fptilt, "\n");
452 fprintf(fprotation, "\n");
453 fprintf(fptheta1, "\n");
454 fprintf(fptheta2, "\n");
455 fprintf(fptheta3, "\n");
456 }
457
458 for (i = 0; i < iCA; i++)
459 {
460 copy_rvec(residuehelixaxis[i], residuehelixaxis_tlast[i]);
461 copy_rvec(residuevector[i], residuevector_tlast[i]);
462 copy_rvec(axis3[i], axis3_tlast[i]);
463 }
464
465 teller++;
466 }
467 while (read_next_x(oenv, status, &t, x, box));
468
469 gmx_rmpbc_done(gpbc);
470
471 gmx_ffclose(fpaxis);
472 gmx_ffclose(fpcenter);
473 xvgrclose(fptilt);
474 xvgrclose(fprotation);
475 gmx_ffclose(fprise);
476 gmx_ffclose(fpradius);
477 gmx_ffclose(fptwist);
478 gmx_ffclose(fpbending);
479 gmx_ffclose(fptheta1);
480 gmx_ffclose(fptheta2);
481 gmx_ffclose(fptheta3);
482
483 close_trx(status);
484
485 return 0;
486 }
The ultimate molecular dynamics simulation package