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