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