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