src/gromacs/gmxana/gmx_helixorient.c

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