Project

General

Profile

gmx_do_dssp.cpp

Boris Timofeev, 08/07/2018 05:17 PM

 
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) 2012,2013,2014,2015,2017 by the GROMACS development team, led by
7
 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8
 * and including many others, as listed in the AUTHORS file in the
9
 * top-level source directory and at http://www.gromacs.org.
10
 *
11
 * GROMACS is free software; you can redistribute it and/or
12
 * modify it under the terms of the GNU Lesser General Public License
13
 * as published by the Free Software Foundation; either version 2.1
14
 * of the License, or (at your option) any later version.
15
 *
16
 * GROMACS is distributed in the hope that it will be useful,
17
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19
 * Lesser General Public License for more details.
20
 *
21
 * You should have received a copy of the GNU Lesser General Public
22
 * License along with GROMACS; if not, see
23
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
25
 *
26
 * If you want to redistribute modifications to GROMACS, please
27
 * consider that scientific software is very special. Version
28
 * control is crucial - bugs must be traceable. We will be happy to
29
 * consider code for inclusion in the official distribution, but
30
 * derived work must not be called official GROMACS. Details are found
31
 * in the README & COPYING files - if they are missing, get the
32
 * official version at http://www.gromacs.org.
33
 *
34
 * To help us fund GROMACS development, we humbly ask that you cite
35
 * the research papers on the package. Check out http://www.gromacs.org.
36
 */
37
#include "gmxpre.h"
38

    
39
#include <cstdlib>
40
#include <cstring>
41

    
42
#include <algorithm>
43

    
44
#include "gromacs/commandline/pargs.h"
45
#include "gromacs/commandline/viewit.h"
46
#include "gromacs/fileio/confio.h"
47
#include "gromacs/fileio/matio.h"
48
#include "gromacs/fileio/pdbio.h"
49
#include "gromacs/fileio/trxio.h"
50
#include "gromacs/fileio/xvgr.h"
51
#include "gromacs/gmxana/gmx_ana.h"
52
#include "gromacs/gmxana/gstat.h"
53
#include "gromacs/pbcutil/rmpbc.h"
54
#include "gromacs/topology/index.h"
55
#include "gromacs/topology/topology.h"
56
#include "gromacs/utility/arraysize.h"
57
#include "gromacs/utility/cstringutil.h"
58
#include "gromacs/utility/dir_separator.h"
59
#include "gromacs/utility/fatalerror.h"
60
#include "gromacs/utility/futil.h"
61
#include "gromacs/utility/smalloc.h"
62
#include "gromacs/utility/strdb.h"
63

    
64
#if GMX_NATIVE_WINDOWS
65
        #define NULL_DEVICE  "nul"
66
        #define popen  _popen
67
    #define pclose _pclose
68
#else
69
        #define NULL_DEVICE  "/dev/null"
70
#endif
71

    
72
static int strip_dssp(FILE *tapeout, int nres,
73
                      gmx_bool bPhobres[], real t,
74
                      real *acc, FILE *fTArea,
75
                      t_matrix *mat, int average_area[],
76
                      const gmx_output_env_t *oenv)
77
{
78
    static gmx_bool bFirst = TRUE;
79
    static char    *ssbuf;
80
    static int      xsize, frame;
81
    char            buf[STRLEN+1];
82
    char            SSTP;
83
    int             i, nr, iacc, nresidues;
84
    int             naccf, naccb; /* Count hydrophobic and hydrophilic residues */
85
    real            iaccf, iaccb;
86
    t_xpmelmt       c;
87

    
88
    /* Skip header */
89
    do
90
    {
91
        fgets2(buf, STRLEN, tapeout);
92
    }
93
    while (std::strstr(buf, "KAPPA") == nullptr);
94
    if (bFirst)
95
    {
96
        /* Since we also have empty lines in the dssp output (temp) file,
97
         * and some line content is saved to the ssbuf variable,
98
         * we need more memory than just nres elements. To be shure,
99
         * we allocate 2*nres-1, since for each chain there is a
100
         * separating line in the temp file. (At most each residue
101
         * could have been defined as a separate chain.) */
102
        snew(ssbuf, 2*nres-1);
103
    }
104

    
105
    iaccb     = iaccf = 0;
106
    nresidues = 0;
107
    naccf     = 0;
108
    naccb     = 0;
109
    for (nr = 0; (fgets2(buf, STRLEN, tapeout) != nullptr); nr++)
110
    {
111
        if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
112
        {
113
            SSTP = '=';     /* Chain separator sign '=' */
114
        }
115
        else
116
        {
117
            SSTP = buf[16] == ' ' ? '~' : buf[16];
118
        }
119
        ssbuf[nr] = SSTP;
120

    
121
        buf[39] = '\0';
122

    
123
        /* Only calculate solvent accessible area if needed */
124
        if ((nullptr != acc) && (buf[13] != '!'))
125
        {
126
            sscanf(&(buf[34]), "%d", &iacc);
127
            acc[nr] = iacc;
128
            /* average_area and bPhobres are counted from 0...nres-1 */
129
            average_area[nresidues] += iacc;
130
            if (bPhobres[nresidues])
131
            {
132
                naccb++;
133
                iaccb += iacc;
134
            }
135
            else
136
            {
137
                naccf++;
138
                iaccf += iacc;
139
            }
140
            /* Keep track of the residue number (do not count chain separator lines '!') */
141
            nresidues++;
142
        }
143
    }
144
    ssbuf[nr] = '\0';
145

    
146
    if (bFirst)
147
    {
148
        if (0 != acc)
149
        {
150
            fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
151
        }
152

    
153
        sprintf(mat->title, "Secondary structure");
154
        mat->legend[0] = 0;
155
        sprintf(mat->label_x, "%s", output_env_get_time_label(oenv));
156
        sprintf(mat->label_y, "Residue");
157
        mat->bDiscrete = TRUE;
158
        mat->ny        = nr;
159
        snew(mat->axis_y, nr);
160
        for (i = 0; i < nr; i++)
161
        {
162
            mat->axis_y[i] = i+1;
163
        }
164
        mat->axis_x = nullptr;
165
        mat->matrix = nullptr;
166
        xsize       = 0;
167
        frame       = 0;
168
        bFirst      = FALSE;
169
    }
170
    if (frame >= xsize)
171
    {
172
        xsize += 10;
173
        srenew(mat->axis_x, xsize);
174
        srenew(mat->matrix, xsize);
175
    }
176
    mat->axis_x[frame] = t;
177
    snew(mat->matrix[frame], nr);
178
    c.c2 = 0;
179
    for (i = 0; i < nr; i++)
180
    {
181
        c.c1                  = ssbuf[i];
182
        mat->matrix[frame][i] = std::max(static_cast<t_matelmt>(0), searchcmap(mat->nmap, mat->map, c));
183
    }
184
    frame++;
185
    mat->nx = frame;
186

    
187
    if (fTArea)
188
    {
189
        fprintf(fTArea, "%10g  %10g  %10g\n", t, 0.01*iaccb, 0.01*iaccf);
190
    }
191

    
192
    /* Return the number of lines found in the dssp file (i.e. number
193
     * of redidues plus chain separator lines).
194
     * This is the number of y elements needed for the area xpm file */
195
    return nr;
196
}
197

    
198
static gmx_bool *bPhobics(t_atoms *atoms)
199
{
200
    int     j,   i, nb;
201
    char      **cb;
202
    gmx_bool   *bb;
203

    
204
        int    n, n_surf;
205
        char    surffn[] = "surface.dat";
206
        char  **surf_res, **surf_lines;
207

    
208
    nb = get_lines("phbres.dat", &cb);
209
    snew(bb, atoms->nres);
210

    
211

    
212
        n_surf = get_lines(surffn, &surf_lines);
213
        snew(surf_res, n_surf);
214
        for (i = 0; (i < n_surf); i++)
215
        {
216
                snew(surf_res[i], 5);
217
                sscanf(surf_lines[i], "%s", surf_res[i]);
218
        }
219

    
220
        for (i = 0, j = 0; (i < atoms->nres); i++)
221
        {
222
                if (-1 != search_str(n_surf, surf_res, *atoms->resinfo[i].name) )
223
                {
224
                        bb[j++] = (-1 != search_str(nb, cb, *atoms->resinfo[i].name)) ? TRUE : FALSE;
225
            }
226
    }
227
        if (i != j)
228
        {
229
                fprintf(stderr, "Not all residues was recognized (%d from %d), the result may be inaccurate!\n", j, i);
230
        }
231
        for (i = 0; (i < n_surf); i++)
232
        {
233
                sfree(surf_res[i]);
234
        }         
235
        sfree(surf_res);
236

    
237
        return bb;
238
}
239

    
240
static void check_oo(t_atoms *atoms)
241
{
242
    char *OOO;
243

    
244
    int   i;
245

    
246
    OOO = gmx_strdup("O");
247

    
248
    for (i = 0; (i < atoms->nr); i++)
249
    {
250
        if (std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
251
        {
252
            *atoms->atomname[i] = OOO;
253
        }
254
        else if (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
255
        {
256
            *atoms->atomname[i] = OOO;
257
        }
258
        else if (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
259
        {
260
            *atoms->atomname[i] = OOO;
261
        }
262
    }
263
}
264

    
265
static void norm_acc(t_atoms *atoms, int nres,
266
                     real av_area[], real norm_av_area[])
267
{
268
    int     i, n, n_surf;
269

    
270
    char    surffn[] = "surface.dat";
271
    char  **surf_res, **surf_lines;
272
    double *surf;
273

    
274
    n_surf = get_lines(surffn, &surf_lines);
275
    snew(surf, n_surf);
276
    snew(surf_res, n_surf);
277
    for (i = 0; (i < n_surf); i++)
278
    {
279
        snew(surf_res[i], 5);
280
        sscanf(surf_lines[i], "%s %lf", surf_res[i], &surf[i]);
281
    }
282

    
283
    for (i = 0; (i < nres); i++)
284
    {
285
        n = search_str(n_surf, surf_res, *atoms->resinfo[i].name);
286
        if (n != -1)
287
        {
288
            norm_av_area[i] = av_area[i] / surf[n];
289
        }
290
        else
291
        {
292
            fprintf(stderr, "Residue %s not found in surface database (%s)\n",
293
                    *atoms->resinfo[i].name, surffn);
294
        }
295
    }
296
}
297

    
298
static void prune_ss_legend(t_matrix *mat)
299
{
300
    gmx_bool  *present;
301
    int       *newnum;
302
    int        i, r, f, newnmap;
303
    t_mapping *newmap;
304

    
305
    snew(present, mat->nmap);
306
    snew(newnum, mat->nmap);
307

    
308
    for (f = 0; f < mat->nx; f++)
309
    {
310
        for (r = 0; r < mat->ny; r++)
311
        {
312
            present[mat->matrix[f][r]] = TRUE;
313
        }
314
    }
315

    
316
    newnmap = 0;
317
    newmap  = nullptr;
318
    for (i = 0; i < mat->nmap; i++)
319
    {
320
        newnum[i] = -1;
321
        if (present[i])
322
        {
323
            newnum[i] = newnmap;
324
            newnmap++;
325
            srenew(newmap, newnmap);
326
            newmap[newnmap-1] = mat->map[i];
327
        }
328
    }
329
    if (newnmap != mat->nmap)
330
    {
331
        mat->nmap = newnmap;
332
        mat->map  = newmap;
333
        for (f = 0; f < mat->nx; f++)
334
        {
335
            for (r = 0; r < mat->ny; r++)
336
            {
337
                mat->matrix[f][r] = newnum[mat->matrix[f][r]];
338
            }
339
        }
340
    }
341
}
342

    
343
static void write_sas_mat(const char *fn, real **accr, int nframe, int nres, t_matrix *mat)
344
{
345
    real  lo, hi;
346
    int   i, j, nlev;
347
    t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
348
    FILE *fp;
349

    
350
    if (fn)
351
    {
352
        hi = lo = accr[0][0];
353
        for (i = 0; i < nframe; i++)
354
        {
355
            for (j = 0; j < nres; j++)
356
            {
357
                lo = std::min(lo, accr[i][j]);
358
                hi = std::max(hi, accr[i][j]);
359
            }
360
        }
361
        fp   = gmx_ffopen(fn, "w");
362
        nlev = static_cast<int>(hi-lo+1);
363
        write_xpm(fp, 0, "Solvent Accessible Surface", "Surface (A^2)",
364
                  "Time", "Residue Index", nframe, nres,
365
                  mat->axis_x, mat->axis_y, accr, lo, hi, rlo, rhi, &nlev);
366
        gmx_ffclose(fp);
367
    }
368
}
369

    
370
static void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
371
                const gmx_output_env_t *oenv)
372
{
373
    FILE        *fp;
374
    t_mapping   *map;
375
    int          f, r, *count, *total, ss_count, total_count;
376
    size_t       s;
377
    const char** leg;
378

    
379
    map = mat->map;
380
    snew(count, mat->nmap);
381
    snew(total, mat->nmap);
382
    snew(leg, mat->nmap+1);
383
    leg[0] = "Structure";
384
    for (s = 0; s < (size_t)mat->nmap; s++)
385
    {
386
        leg[s+1] = gmx_strdup(map[s].desc);
387
    }
388

    
389
    fp = xvgropen(outfile, "Secondary Structure",
390
                  output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);
391
    if (output_env_get_print_xvgr_codes(oenv))
392
    {
393
        fprintf(fp, "@ subtitle \"Structure = ");
394
    }
395
    for (s = 0; s < std::strlen(ss_string); s++)
396
    {
397
        if (s > 0)
398
        {
399
            fprintf(fp, " + ");
400
        }
401
        for (f = 0; f < mat->nmap; f++)
402
        {
403
            if (ss_string[s] == map[f].code.c1)
404
            {
405
                fprintf(fp, "%s", map[f].desc);
406
            }
407
        }
408
    }
409
    fprintf(fp, "\"\n");
410
    xvgr_legend(fp, mat->nmap+1, leg, oenv);
411

    
412
    total_count = 0;
413
    for (s = 0; s < (size_t)mat->nmap; s++)
414
    {
415
        total[s] = 0;
416
    }
417
    for (f = 0; f < mat->nx; f++)
418
    {
419
        ss_count = 0;
420
        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
421
        {
422
            count[s] = 0;
423
        }
424
        for (r = 0; r < mat->ny; r++)
425
        {
426
            count[mat->matrix[f][r]]++;
427
            total[mat->matrix[f][r]]++;
428
        }
429
        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
430
        {
431
            if (std::strchr(ss_string, map[s].code.c1))
432
            {
433
                ss_count    += count[s];
434
                total_count += count[s];
435
            }
436
        }
437
        fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
438
        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
439
        {
440
            fprintf(fp, " %5d", count[s]);
441
        }
442
        fprintf(fp, "\n");
443
    }
444
    /* now print column totals */
445
    fprintf(fp, "%-8s %5d", "# Totals", total_count);
446
    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
447
    {
448
        fprintf(fp, " %5d", total[s]);
449
    }
450
    fprintf(fp, "\n");
451

    
452
    /* now print percentages */
453
    fprintf(fp, "%-8s %5.2f", "# SS %", total_count / static_cast<real>(mat->nx * mat->ny));
454
    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
455
    {
456
        fprintf(fp, " %5.2f", total[s] / static_cast<real>(mat->nx * mat->ny));
457
    }
458
    fprintf(fp, "\n");
459

    
460
    xvgrclose(fp);
461
    sfree(leg);
462
    sfree(count);
463
}
464

    
465
int gmx_do_dssp(int argc, char *argv[])
466
{
467
    const char        *desc[] = {
468
        "[THISMODULE] ",
469
        "reads a trajectory file and computes the secondary structure for",
470
        "each time frame ",
471
        "calling the dssp program. If you do not have the dssp program,",
472
        "get it from http://swift.cmbi.ru.nl/gv/dssp. [THISMODULE] assumes ",
473
        "that the dssp executable is located in ",
474
        "[TT]/usr/local/bin/dssp[tt]. If this is not the case, then you should",
475
        "set an environment variable [TT]DSSP[tt] pointing to the dssp",
476
        "executable, e.g.: [PAR]",
477
        "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
478
        "Since version 2.0.0, dssp is invoked with a syntax that differs",
479
        "from earlier versions. If you have an older version of dssp,",
480
        "use the [TT]-ver[tt] option to direct do_dssp to use the older syntax.",
481
        "By default, do_dssp uses the syntax introduced with version 2.0.0.",
482
        "Even newer versions (which at the time of writing are not yet released)",
483
        "are assumed to have the same syntax as 2.0.0.[PAR]",
484
        "The structure assignment for each residue and time is written to an",
485
        "[REF].xpm[ref] matrix file. This file can be visualized with for instance",
486
        "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
487
        "Individual chains are separated by light grey lines in the [REF].xpm[ref] and",
488
        "postscript files.",
489
        "The number of residues with each secondary structure type and the",
490
        "total secondary structure ([TT]-sss[tt]) count as a function of",
491
        "time are also written to file ([TT]-sc[tt]).[PAR]",
492
        "Solvent accessible surface (SAS) per residue can be calculated, both in",
493
        "absolute values (A^2) and in fractions of the maximal accessible",
494
        "surface of a residue. The maximal accessible surface is defined as",
495
        "the accessible surface of a residue in a chain of glycines.",
496
        "[BB]Note[bb] that the program [gmx-sas] can also compute SAS",
497
        "and that is more efficient.[PAR]",
498
        "Finally, this program can dump the secondary structure in a special file",
499
        "[TT]ssdump.dat[tt] for usage in the program [gmx-chi]. Together",
500
        "these two programs can be used to analyze dihedral properties as a",
501
        "function of secondary structure type."
502
    };
503
    static gmx_bool    bVerbose;
504
    static const char *ss_string   = "HEBT";
505
    static int         dsspVersion = 2;
506
    t_pargs            pa[]        = {
507
        { "-v",  FALSE, etBOOL, {&bVerbose},
508
          "HIDDENGenerate miles of useless information" },
509
        { "-sss", FALSE, etSTR, {&ss_string},
510
          "Secondary structures for structure count"},
511
        { "-ver", FALSE, etINT, {&dsspVersion},
512
          "DSSP major version. Syntax changed with version 2"}
513
    };
514

    
515
    t_trxstatus       *status;
516
    FILE              *tapein, *tapeout;
517
    FILE              *ss, *acc, *fTArea, *tmpf;
518
    const char        *fnSCount, *fnArea, *fnTArea, *fnAArea;
519
    const char        *leg[] = { "Phobic", "Phylic" };
520
    t_topology         top;
521
    int                ePBC;
522
    t_atoms           *atoms;
523
    t_matrix           mat;
524
    int                nres, nr0, naccr, nres_plus_separators;
525
    gmx_bool          *bPhbres, bDoAccSurf;
526
    real               t;
527
    int                i, j, natoms, nframe = 0;
528
    matrix             box = {{0}};
529
    int                gnx;
530
    char              *grpnm, *ss_str;
531
    int               *index;
532
    rvec              *xp, *x;
533
    int               *average_area;
534
    real             **accr, *accr_ptr = nullptr, *av_area, *norm_av_area;
535
    char               pdbfile[32], tmpfile[32];
536
    char               dssp[256];
537
    const char        *dptr;
538
    gmx_output_env_t  *oenv;
539
    gmx_rmpbc_t        gpbc = nullptr;
540

    
541
    t_filenm           fnm[] = {
542
        { efTRX, "-f",   nullptr,      ffREAD },
543
        { efTPS, nullptr,   nullptr,      ffREAD },
544
        { efNDX, nullptr,   nullptr,      ffOPTRD },
545
        { efDAT, "-ssdump", "ssdump", ffOPTWR },
546
        { efMAP, "-map", "ss",      ffLIBRD },
547
        { efXPM, "-o",   "ss",      ffWRITE },
548
        { efXVG, "-sc",  "scount",  ffWRITE },
549
        { efXPM, "-a",   "area",    ffOPTWR },
550
        { efXVG, "-ta",  "totarea", ffOPTWR },
551
        { efXVG, "-aa",  "averarea", ffOPTWR }
552
    };
553
#define NFILE asize(fnm)
554

    
555
    if (!parse_common_args(&argc, argv,
556
                           PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT,
557
                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
558
    {
559
        return 0;
560
    }
561
    fnSCount   = opt2fn("-sc", NFILE, fnm);
562
    fnArea     = opt2fn_null("-a", NFILE, fnm);
563
    fnTArea    = opt2fn_null("-ta", NFILE, fnm);
564
    fnAArea    = opt2fn_null("-aa", NFILE, fnm);
565
    bDoAccSurf = (fnArea || fnTArea || fnAArea);
566

    
567
    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xp, nullptr, box, FALSE);
568
    atoms = &(top.atoms);
569
    check_oo(atoms);
570
    bPhbres = bPhobics(atoms);
571

    
572
    get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpnm);
573
    nres = 0;
574
    nr0  = -1;
575
    for (i = 0; (i < gnx); i++)
576
    {
577
        if (atoms->atom[index[i]].resind != nr0)
578
        {
579
            nr0 = atoms->atom[index[i]].resind;
580
            nres++;
581
        }
582
    }
583
    fprintf(stderr, "There are %d residues in your selected group\n", nres);
584

    
585
    std::strcpy(pdbfile, "ddXXXXXX");
586
    gmx_tmpnam(pdbfile);
587
    if ((tmpf = fopen(pdbfile, "w")) == nullptr)
588
    {
589
        sprintf(pdbfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
590
        gmx_tmpnam(pdbfile);
591
        if ((tmpf = fopen(pdbfile, "w")) == nullptr)
592
        {
593
            gmx_fatal(FARGS, "Can not open tmp file %s", pdbfile);
594
        }
595
    }
596
    else
597
    {
598
        fclose(tmpf);
599
    }
600
        remove(pdbfile); // To exclude empty backup
601

    
602
    std::strcpy(tmpfile, "ddXXXXXX");
603
    gmx_tmpnam(tmpfile);
604
    if ((tmpf = fopen(tmpfile, "w")) == nullptr)
605
    {
606
        sprintf(tmpfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
607
        gmx_tmpnam(tmpfile);
608
        if ((tmpf = fopen(tmpfile, "w")) == nullptr)
609
        {
610
            gmx_fatal(FARGS, "Can not open tmp file %s", tmpfile);
611
        }
612
    }
613
    else
614
    {
615
        fclose(tmpf);
616
    }
617
        remove(tmpfile);  // To exclude empty backup
618
    
619
        if ((dptr = getenv("DSSP")) == nullptr)
620
    {
621
        dptr = "/usr/local/bin/dssp";
622
    }
623
    if (!gmx_fexist(dptr))
624
    {
625
        gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)",
626
                  dptr);
627
    }
628
    if (dsspVersion >= 2)
629
    {
630
        if (dsspVersion > 2)
631
        {
632
            printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by do_dssp. Assuming version 2 syntax.\n\n", dsspVersion);
633
        }
634
#if HAVE_PIPES || GMX_NATIVE_WINDOWS
635
                sprintf(dssp, "%s -i %s %s",
636
                        dptr, pdbfile, bVerbose ? "" : "2>" NULL_DEVICE);
637
#else
638
                sprintf(dssp, "%s -i %s -o %s > %s %s",
639
                        dptr, pdbfile, tmpfile, NULL_DEVICE, bVerbose ? "" : "2>" NULL_DEVICE);
640
#endif
641
    }
642
    else
643
    {
644
#if HAVE_PIPES || GMX_NATIVE_WINDOWS
645
                sprintf(dssp, "%s %s %s %s",
646
                        dptr, bDoAccSurf ? "" : "-na", pdbfile, bVerbose ? "" : "2>" NULL_DEVICE);
647
#else
648
                sprintf(dssp, "%s %s %s %s > %s %s",
649
                        dptr, bDoAccSurf ? "" : "-na", pdbfile, tmpfile, NULL_DEVICE, bVerbose ? "" : "2>" NULL_DEVICE);
650
#endif
651
        }
652
    fprintf(stderr, "dssp cmd='%s'\n", dssp);
653

    
654
    if (fnTArea)
655
    {
656
        fTArea = xvgropen(fnTArea, "Solvent Accessible Surface Area",
657
                          output_env_get_xvgr_tlabel(oenv), "Area (nm\\S2\\N)", oenv);
658
        xvgr_legend(fTArea, 2, leg, oenv);
659
    }
660
    else
661
    {
662
        fTArea = nullptr;
663
    }
664

    
665
    mat.map  = nullptr;
666
    mat.nmap = readcmap(opt2fn("-map", NFILE, fnm), &(mat.map));
667

    
668
    natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
669
    if (natoms > atoms->nr)
670
    {
671
        gmx_fatal(FARGS, "\nTrajectory does not match topology (%d > %d)!", natoms, atoms->nr);
672
    }
673
    if (gnx > natoms)
674
    {
675
        gmx_fatal(FARGS, "\nTrajectory does not match selected group!");
676
    }
677

    
678
    snew(average_area, atoms->nres);
679
    snew(av_area, atoms->nres);
680
    snew(norm_av_area, atoms->nres);
681
    accr  = nullptr;
682
    naccr = 0;
683

    
684
    gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
685
    do
686
    {
687
        t = output_env_conv_time(oenv, t);
688
        if (bDoAccSurf && nframe >= naccr)
689
        {
690
            naccr += 10;
691
            srenew(accr, naccr);
692
            for (i = naccr-10; i < naccr; i++)
693
            {
694
                snew(accr[i], 2*atoms->nres-1);
695
            }
696
        }
697
        gmx_rmpbc(gpbc, natoms, box, x);
698
        tapein = gmx_ffopen(pdbfile, "w");
699
        write_pdbfile_indexed(tapein, nullptr, atoms, x, ePBC, box, ' ', -1, gnx, index, nullptr, TRUE);
700
        gmx_ffclose(tapein);
701
#if !HAVE_PIPES && !GMX_NATIVE_WINDOWS
702
        if (0 != system(dssp))
703
        {
704
            gmx_fatal(FARGS, "Failed to execute command: %s\n",
705
                      "Try specifying your dssp version with the -ver option.", dssp);
706
        }
707
#endif
708

    
709
        /* strip_dssp returns the number of lines found in the dssp file, i.e.
710
         * the number of residues plus the separator lines */
711

    
712
        if (bDoAccSurf)
713
        {
714
            accr_ptr = accr[nframe];
715
        }
716
#if HAVE_PIPES || GMX_NATIVE_WINDOWS
717
                if ( nullptr == (tapeout = popen(dssp, "r")))
718
                {
719
                        gmx_fatal(FARGS, "Failed to execute command: %s\n",
720
                                "Try specifying your dssp version with the -ver option.", dssp);
721
                }
722
#else
723
                tapeout = gmx_ffopen(tmpfile, "r");
724
#endif
725
        nres_plus_separators = strip_dssp(tapeout, nres, bPhbres, t,
726
                                          accr_ptr, fTArea, &mat, average_area, oenv);
727
#if HAVE_PIPES || GMX_NATIVE_WINDOWS
728
                pclose(tapeout);
729
#else
730
                gmx_ffclose(tapeout);
731
        remove(tmpfile);
732
#endif
733
        remove(pdbfile);
734
        nframe++;
735
    }
736
    while (read_next_x(oenv, status, &t, x, box));
737

    
738
    fprintf(stderr, "\n");
739
    close_trj(status);
740
    if (fTArea)
741
    {
742
        xvgrclose(fTArea);
743
    }
744
    gmx_rmpbc_done(gpbc);
745

    
746
    prune_ss_legend(&mat);
747

    
748
    ss        = opt2FILE("-o", NFILE, fnm, "w");
749
    mat.flags = 0;
750
    write_xpm_m(ss, mat);
751
    gmx_ffclose(ss);
752

    
753
    if (opt2bSet("-ssdump", NFILE, fnm))
754
    {
755
        ss = opt2FILE("-ssdump", NFILE, fnm, "w");
756
        snew(ss_str, nres+1);
757
        fprintf(ss, "%d\n", nres);
758
        for (j = 0; j < mat.nx; j++)
759
        {
760
            for (i = 0; (i < mat.ny); i++)
761
            {
762
                ss_str[i] = mat.map[mat.matrix[j][i]].code.c1;
763
            }
764
            ss_str[i] = '\0';
765
            fprintf(ss, "%s\n", ss_str);
766
        }
767
        gmx_ffclose(ss);
768
        sfree(ss_str);
769
    }
770
    analyse_ss(fnSCount, &mat, ss_string, oenv);
771

    
772
    if (bDoAccSurf)
773
    {
774
        write_sas_mat(fnArea, accr, nframe, nres_plus_separators, &mat);
775

    
776
        for (i = 0; i < atoms->nres; i++)
777
        {
778
            av_area[i] = (average_area[i] / static_cast<real>(nframe));
779
        }
780

    
781
        norm_acc(atoms, nres, av_area, norm_av_area);
782

    
783
        if (fnAArea)
784
        {
785
            acc = xvgropen(fnAArea, "Average Accessible Area",
786
                           "Residue", "A\\S2", oenv);
787
            for (i = 0; (i < nres); i++)
788
            {
789
                fprintf(acc, "%5d  %10g %10g\n", i+1, av_area[i], norm_av_area[i]);
790
            }
791
            xvgrclose(acc);
792
        }
793
    }
794

    
795
    view_all(oenv, NFILE, fnm);
796

    
797
    return 0;
798
}