src/contrib/dsspcore.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_dsspcore_c = "$Id$";
  30
  31 /* Output from p2c, the Pascal-to-C translator */
  32 /* From input file "dssp.p" */
  33
  34
  35
  36 /* ------------------------------------------------------------------ */
  37 /*
  38
  39 DSSP version October 1988.
  40 This copy for Rudi Drunen at Univ_Groningen
  41 who have agreed to the following software license agreement:
  42
  43
  44 An academic license for the DSSP program
  45 ((c) W. Kabsch, C. Sander and MPI-MF, 1983, 1985, 1988)
  46 is granted to in exchange for the following commitments:
  47
  48 I hereby certify that
  49
  50         (1) I am an academic user at an academic research institution. In
  51             using the software, I will respect the interests of the authors
  52             and their institutions.
  53
  54         (2) I will not use the software in commercial activities without
  55             a written commercial license agreement; commercial activities
  56             include, in particular, work under contract from a commercial
  57             company.
  58
  59         (3) I will not redistribute the software to others outside of my
  60             immediate research group. I will suggest to other interested
  61             research groups to contact the authors directly.
  62
  63         (4) I will not alter or suppress the run-time copyright message.
  64
  65         (5) I will acknowledge the program authors on any publication of
  66             scientific results based in part on use of the program and
  67             cite the article in which the program was described.
  68
  69         (6) I will report evidence of program bugs to the authors.
  70
  71         (7) I will send the source code of any bug corrections and program
  72             extensions, major or minor, to the original authors, for free
  73             academic use. If I have made major extensions which are incor-
  74             porated by the authors, I reserve the right to be appropriately
  75             included in any future commercial license agreement.
  76
  77         (8) I will not extract part of the software, e.g. modules or sub-
  78             routines, for use in other contexts without permission by the
  79             authors.
  80
  81         (9) I will not use the program in the context of classified research.
  82 */
  83 /* PREVIOUS RELEASE: VERSION OCTOBER 1985                             */
  84 /* PREVIOUS RELEASE: VERSION JUNE 1983                                */
  85 /* LANGUAGE: STANDARD PASCAL WITH 128 CHARACTER ASCII SET             */
  86 /* AUTHORS AND COPYRIGHT (1983,1985,1988):
  87    Wolfgang Kabsch and Chris Sander, Max Planck Institut
  88    fuer Medizinische Forschung, Jahnstr. 29, 6900 Heidelberg, Germany
  89    Telephone: +49-6221-486 276  Telex: 461505 mpimf d
  90    Bitnet:    KABSCH@EMBL
  91    Current address for Chris Sander:
  92    Biocomputing, EMBL, 6900 Heidelberg, Germany
  93    Telephone: +49-6221-387 361 Telex: 461613 embl d
  94    Telefax:   +49-6221-387 306
  95    Bitnet:    SANDER@EMBL
  96    Do report errors if you find any.
  97    Reference: Kabsch,W. and Sander,C. (1983) Biopolymers 22, 2577-2637*/
  98 /*--------------------------------------------------------------------*/
  99 /* DEFINES SECONDARY STRUCTURE AND SOLVENT EXPOSURE OF PROTEINS FROM
 100    ATOMIC COORDINATES AS GIVEN BY THE BROOKHAVEN PROTEIN DATA BANK.   */
 101 /*--------------------------------------------------------------------*/
 102 /* This program including sample input and output files for dataset 1PPT
 103    is available from the authors in exchange for an academic or
 104    commercial license agreement. The program is no longer available
 105    from the Brookhaven Protein Data Bank */
 106 /*--------------------------------------------------------------------*/
 107 /* CORRECTION AND MODIFICATION LOG SINCE JUNE 1983 */
 108 /* (1) MODIFICATIONS THAT AFFECT OUTPUT ON FILE TAPEOUT FOR AT LEAST ONE
 109        OF THE 62 PROTEIN DATA SETS IN THE 1983 BIOPOLYMERS PAPER:
 110    - SIDECHAIN ATOMS MORE THAN MAXDIST ANGSTROM DISTANT FROM ATOM CA ARE
 111      DECLARED ILLEGAL AND IGNORED. OUTPUT CHANGE: ACCESSIBILITY VALUES
 112      FOR ASN 76 OF 1SBT (ILLEGAL ATOM OD1) AND PRO 49 OF 156B (ILLEGAL
 113      ATOM UNK).
 114    - ANY RESIDUE WITH INCOMPLETE BACKBONE IS IGNORED. OUTPUT CHANGE:
 115      CHAIN BREAK BETWEEN RESIDUE SER 11 AND ILE 16 IN 2GCH
 116      (DUE TO INCOMPLETE COORDINATES FOR SER 11) IS NOW CORRECT.
 117    (2) MODIFICATIONS THAT DO NOT AFFECT OUTPUT ON FILE TAPEOUT FOR ANY
 118        OF THE 62 PROTEIN DATA SETS IN THE 1983 BIOPOLYMERS PAPER:
 119    - SPELLING OF FLAGCHIRALITY AND TESTSSBOND CORRECTED.
 120    - WARNING MESSAGE FOR RESIDUES WITH NON-STANDARD NUMBER OF
 121      SIDECHAIN ATOMS HAS BEEN ADDED. FOR EXAMPLE, THIS ALERTS THE USER
 122      TO BAD DATA FOR RESIDUES 8,12,21,24 AND 44 OF DATA SET 2RXN.
 123    - WARNING MESSAGE FOR RESIDUES IGNORED DUE TO NON-STANDARD RESIDUE
 124      NAME SUCH AS 'ACE' AND 'FOR' HAS BEEN ADDED.
 125    - WARNING MESSAGE FOR ALTERNATE ATOM LOCATION IDENTIFIER HAS BEEN
 126      ADDED. FOR EXAMPLE, THE USER IS NOW WARNED THAT ATOM CH2 IN ANY
 127      TRP OF DATA SET 1APP IS IGNORED DUE TO BAD ALTERNATE LOCATION
 128      IDENTIFIER '2'.
 129    WE THANK STEVEN SHERIFF, FRANCOIS COLONNA AND JOHN MOULT FOR
 130    REPORTING PROBLEMS AND STEVEN SHERIF, JANET THORNTON AND
 131    WILLIE TAYLOR FOR RESULTS OF TEST RUNS ON VAX COMPUTERS.
 132
 133    Changes after 1985:
 134
 135    - program speeded up by a factor of two or three by avoiding use
 136      of square root.
 137    - hydrogen atoms in data set ignored on input (H of NH of backbone
 138      is built as before)
 139    - 18-AUG-1988: CADIST=9.0, replacing CADIST=8.0. Has affected output
 140      for 63/300 proteins in a minor way. Thanks to Jean Richelle (Bruxelles)
 141      for pointing out this bug.
 142
 143      Output changes due to change in parameter CADIST (8 to 9 Angstrom) :
 144      additional backbone-backbone Hbonds found with slight
 145      adjustments in secondary structure summary. In about 300 protein
 146      data sets from the Fall 1988 PDB release, 63 additional
 147      Hbonds were found, i.e. 0.2 Hbonds per protein (29 of type
 148      i,i+5;  16 of type i,i+4; 6 of type i,i+3; 10 in antiparallel beta
 149      bridges and 2 in a parallel beta bridge). These additional
 150      Hbonds affected the secondary structure summary of 26 of these
 151      protein data sets in a minor way, typically joining a 4-turn to
 152      an alpha-helix, changing a geometrical turn to a hydrogen-
 153      bonded turn or adding an extra residue pair to a beta ladder.
 154      The changes are (using _ for blank):
 155
 156      [protein id, old secstruc summary > corrected summary]
 157
 158      1FC2_E > EE   and  _T > ET
 159      1GP1GGG > HHH
 160      1HBSS > T
 161      1HDSS > T and  GGGGGG > TTHHHH
 162      1HFM__ > TT
 163      1HKGSSS > TTT
 164      1IG2S_ > TT
 165      1LDX GGG > HTT
 166      1MEV__ > TT  and  _BS > TBS  and  SSS > TTS
 167      1PFCSSS > TTS
 168      1PP2_E > EE  and  _S > ES
 169      1RN3E_SS_E > EEEEEE  and _E > EE  (>3-res beta bulge)
 170      1RNSsame as 1RN3
 171      2ATCHH > TT
 172      2CABB_ > EE
 173      2CPPSS > TT  and  GGGGGG > HHHHTT
 174      2LYZT > H
 175      2MDHSSS > TTT
 176      3CPA TTT > HHH
 177      4CATTTT > HHH
 178      4SBVS > T
 179      5API_ > B
 180      5CPATTT > HHH
 181      7LYZS > H
 182      8CAT_ > B  and  _ > B
 183      8LYZT > H
 184
 185      Note that this bugfix results in a small variation in the total
 186      number of Hbonds, compared to the variation which would
 187      result, say, from changing the (somewhat arbitrary) cutoff of
 188      -0.5 kcal/mol for the Hbond electrostatic potential energy. We
 189      cannot here solve the fundamental difficulty of arbitrary
 190      cutoffs involved in extracting binary classifications (an Hbond
 191      exists, yes/no) from real numbers (coordinates, energies).
 192      However, for most purposes the secondary structure summary agrees
 193      will with anyone's intuitive definition, especially for well-refined and
 194      high resolution structures. For a more clearcut assignment of protein
 195      substructure, we recommend using the detailed H-bond and other assignments
 196      in the columns following the summary column, i.e. columns 19-38 (xxx):
 197
 198      ....;....1....;....2....;....3....;....4....;....5....;....6....;....7..
 199                        xxxxxxxxxxxxxxxxxxxx
 200                        .-- 3-turns/helix
 201                        |.-- 4-turns/helix
 202                        ||.-- 5-turns/helix
 203                        |||.-- geometrical bend
 204                        ||||.-- chirality
 205                        |||||.-- beta bridge label
 206                        ||||||.-- beta bridge label
 207                        |||||||   .-- beta bridge partner resnum
 208                        |||||||   |   .-- beta bridge partner resnum
 209                        |||||||   |   |.-- beta sheet label
 210                        |||||||   |   ||   .-- solvent accessibility
 211                        |||||||   |   ||   |
 212         35   47   I  E     +     0   0    2
 213         36   48   R  E >  S- K   0  39C  97
 214         37   49   Q  T 3  S+     0   0   86    (example from 1EST)
 215         38   50   N  T 3  S+     0   0   34
 216         39   51   W  E <   -KL  36  98C   6
 217                                                                            */
 218 /*--------------------------------------------------------------------*/
 219 /* GENERAL PROGRAM INSTALLATION GUIDE. */
 220 /* (1) THE PROGRAM REQUIRES THE FULL STANDARD ASCII 128 CHARACTER SET,
 221        IN PARTICULAR LOWER CASE LETTERS 'abcdefg....'.
 222    (2) STANDARD PASCAL MAY NOT RECOGNIZE REAL NUMBERS SUCH AS .1, +.1,
 223        -.1 ON INPUT. CHANGE TO 0.1,+0.1,-0.1.
 224    (3) THE NON-STANDARD PROCEDURE 'DATE' RETURNS THE CURRENT DAY, MONTH,
 225        AND YEAR. IF THE PROCEDURE IS NOT CALLED (LINE COMMENTED OUT)
 226        THE PSYCHEDELIC DATE DEC 24, 2001 IS RETURNED. YOU MAY  REPLACE
 227        'DATE' BY THE CORRESPONDING PROCEDURE FROM YOUR PASCAL
 228        IMPLEMENTATION. THE EXAMPLE GIVEN WORKS IN DEC VAX VMS 5.0.
 229    (4) DUE TO INCOMPATIBLE ASCII CODES, SQUARE BRACKETS '[' AND ']'
 230        MAY APPEAR AS '!','?' ETC. USE YOUR EDITOR TO CONVERT THESE.   */
 231 /* INSTALLATION GUIDE FOR VAX/VMS USERS. */
 232 /* (1) THE /OPTIMIZE OPTION OF THE PASCAL COMPILER PRODUCED
 233        INCORRECT CODE ON THE VAX 8600 AT EMBL RUNNING UNDER VMS V4.2.
 234        LATER VERSIONS OF VMS (E.G. VMS 5.0) PRODUCED CORRECT CODE.
 235        IF IN DOUBT, COMPILE USING PASCAL /NOOPTIMIZE.
 236    (2) COPY BROOKHAVEN DATA BANK COORDINATE INPUT TO A FILE NAMED
 237        TAPEIN.DAT . OUTPUT WILL BE IN A FILE NAMED TAPEOUT.DAT        */
 238 /* IMPLEMENTATION ON OTHER COMPUTERS */
 239 /* (1) NORD-500. EXECUTION TIME COMPARABLE TO VAX 780.
 240    (2) SUN-3.    EXECUTION TIME COMPARABLE TO VAX 780.
 241                  Compile using: pc -L
 242                  in ORDER to map upper case letters in keywords
 243                  and identifiers to lower case.
 244    (3) ATARI 520 ST. RUNS FACTOR 60 SLOWER THAN NORD-500 DUE TO
 245        SOFTWARE-EMULATED FLOATING POINT OPERATIONS ON MC68000.        */
 246 /*--------------------------------------------------------------------*/
 247 /* INPUT/OUTPUT FILES. */
 248 /* INPUT:   DEFAULT  INPUT UNIT, E.G. YOUR TERMINAL
 249    OUTPUT:  DEFAULT OUTPUT UNIT, E.G. YOUR TERMINAL,
 250             USED FOR RUN-TIME MESSAGES. WARNINGS AND ERRORS LOOK
 251             LIKE THIS: !!! TEXT !!!
 252    TAPEIN:  FILE WITH PROTEIN DATA BANK COORDINATES, E.G. PDB3PTI.COO
 253    TAPEOUT: DSSP OUTPUT OF LINE LENGTH 128, E.G. PAPER PRINTER        */
 254 /*--------------------------------------------------------------------*/
 255 /* DESCRIPTION OF OUTPUT ON FILE TAPEOUT:
 256    LINE LENGTH OF OUTPUT IS 128 CHARCTERS.
 257    FOR DEFINITONS, SEE ABOVE BIOPOLYMERS ARTICLE.
 258    IN ADDITION NOTE:
 259    HISTOGRAMS - E.G. 2 UNDER COLUMN '8' IN LINE 'RESIDUES PER ALPHA
 260             HELIX' MEANS: THERE ARE 2 ALPHA HELICES OF LENGTH  8
 261             RESIDUES IN THIS DATA SET.
 262    #  RESIDUE AA STRUCTURE BP1 BP2 ACC ..ETC..FOR EACH RESIDUE I:
 263    #  RESIDUE - TWO COLUMNS OF RESIDUE NUMBERS. FIRST COLUMN IS DSSP'S
 264             SEQUENTIAL RESIDUE NUMBER, STARTING AT THE FIRST
 265             RESIDUE ACTUALLY IN THE DATA SET AND INCLUDING CHAIN BREAKS;
 266             THIS NUMBER IS USED TO REFER TO RESIDUES THROUGHOUT. SECOND
 267             COLUMN GIVES CRYSTALLOGRAPHERS' 'RESIDUE SEQUENCE
 268             NUMBER','INSERTION CODE' AND 'CHAIN IDENTIFIER' (SEE PROTEIN
 269             DATA BANK FILE RECORD FORMAT MANUAL), GIVEN FOR REFERENCE
 270             ONLY AND NOT USED FURTHER..
 271    AA -     ONE LETTER AMINO ACID CODE, LOWER CASE FOR SS-BRIDGE CYS.
 272    STRUCTURE - SEE BIOPOLYMERS
 273    BP1 BP2  - RESIDUE NUMBER OF FIRST AND SECOND BRIDGE PARTNER
 274             FOLLOWED BY ONE LETTER SHEET LABEL
 275    ACC -    NUMBER OF WATER MOLECULES IN CONTACT WITH THIS RESIDUE *10.
 276             OR RESIDUE WATER EXPOSED SURFACE IN ANGSTROM**2.
 277    N-H-->O ETC. -  HYDROGEN BONDS. E.G. -3,-1.4 MEANS: IF THIS RESIDUE
 278             IS RESIDUE I THEN N-H OF I IS H-BONDED TO C=O OF I-3
 279             WITH AN ELECTROSTATIC H-BOND ENERGY OF -1.4 KCAL/MOL.
 280    TCO -    COSINE OF ANGLE BETWEEN C=O OF RESIDUE I AND C=O OF
 281             RESIDUE I-1. FOR ALPHA-HELICES, TCO IS NEAR +1, FOR
 282             BETA-SHEETS TCO IS NEAR -1. NOT USED FOR STRUCTURE
 283             DEFINITION.
 284    KAPPA -  VIRTUAL BOND ANGLE (BEND ANGLE) DEFINED BY THE THREE
 285             C-ALPHA ATOMS OF RESIDUES I-2,I,I+2. USED TO DEFINE
 286             BEND (STRUCTURE CODE 'S').
 287    ALPHA -  VIRTUAL TORSION ANGLE (DIHEDRAL ANGLE) DEFINED BY THE FOUR
 288             C-ALPHA ATOMS OF RESIDUES I-1,I,I+1,I+2. USED TO DEFINE
 289             CHIRALITY (STRUCTURE CODE '+' OR '-').
 290    PHI PSI - IUPAC PEPTIDE BACKBONE TORSION ANGLES
 291    X-CA Y-CA Z-CA -  ECHO OF C-ALPHA ATOM COORDINATES              */
 292 /*--------------------------------------------------------------------*/
 293 /* WORDS OF CAUTION */
 294 /* THE VALUES FOR SOLVENT EXPOSURE MAY NOT MEAN WHAT YOU THINK!
 295     (A) EFFECTS LEADING TO LARGER THAN EXPECTED VALUES:
 296      SOLVENT EXPOSURE CALCULATION IGNORES UNUSUAL RESIDUES, LIKE ACE,
 297      OR RESIDUES WITH INCOMPLETE BACKBONE, LIKE ALA 1 OF DATA SET 1CPA.
 298      IT ALSO IGNORES HETATOMS, LIKE A HEME OR METAL LIGANDS.
 299      ALSO, SIDE CHAINS MAY BE INCOMPLETE (AN ERROR MESSAGE IS WRITTEN).
 300     (B) EFFECTS LEADING TO SMALLER THAN EXPECTED VALUES:
 301      IF YOU APPLY THIS PROGRAM TO PROTEIN DATA BANK DATA SETS
 302      CONTAINING OLIGOMERS, SOLVENT EXPOSURE IS FOR THE ENTIRE ASSEMBLY,
 303      NOT FOR THE MONOMER. ALSO, ATOM OXT OF C-TERMINAL RESIDUES IS
 304      TREATED LIKE A SIDE CHAIN ATOM IF IT IS LISTED AS PART OF THE LAST
 305      RESIDUE. ALSO, PEPTIDE SUBSTRATES, WHEN LISTED AS ATOMS RATHER THAN
 306      HETATOMS, ARE TREATED AS PART OF THE PROTEIN, E.G. RESIDUES 499 S
 307      AND 500 S IN 1CPA.                                               */
 308 /* UNKNOWN OR UNUSUAL RESIDUES ARE NAMED X ON OUTPUT AND THEY ARE
 309    NOT CHECKED FOR STANDARD NUMBER OF SIDECHAIN ATOMS.                */
 310 /* ALL EXPLICIT WATER MOLECULES, LIKE OTHER HETATOMS, ARE IGNORED.    */
 311 /* END OF INTRODUCTORY COMMENTS */
 312 /**********************************************************************/
 313
 314
 315 /************************************************************
 316 *                                                           *
 317 *          p2c.h                                            *
 318 *                                                           *
 319 ************************************************************/
 320
 321
 322 #ifndef P2C_H
 323 #define P2C_H
 324
 325
 326 /* Header file for code generated by "p2c", the Pascal-to-C translator */
 327
 328 /* "p2c"  Copyright (C) 1989, 1990, 1991 Free Software Foundation.
 329  * By Dave Gillespie, daveg@csvax.cs.caltech.edu.  Version 1.19.
 330  * This file may be copied, modified, etc. in any way.  It is not restricted
 331  * by the licence agreement accompanying p2c itself.
 332  */
 333
 334
 335 #include <stdio.h>
 336
 337
 338
 339 /* If the following heuristic fails, compile -DBSD=0 for non-BSD systems,
 340    or -DBSD=1 for BSD systems. */
 341
 342 #ifdef M_XENIX
 343 # define BSD 0
 344 #endif
 345
 346 #ifdef vms
 347 # define BSD 0
 348 # ifndef __STDC__
 349 #  define __STDC__ 1
 350 # endif
 351 #endif
 352
 353 #ifdef __TURBOC__
 354 # define MSDOS 1
 355 #endif
 356
 357 #ifdef MSDOS
 358 # define BSD 0
 359 #endif
 360
 361 #ifdef FILE       /* a #define in BSD, a typedef in SYSV (hp-ux, at least) */
 362 # ifndef BSD      /*  (a convenient, but horrible kludge!) */
 363 #  define BSD 1
 364 # endif
 365 #endif
 366
 367 #ifdef BSD
 368 # if !BSD
 369 #  undef BSD
 370 # endif
 371 #endif
 372
 373
 374 #ifdef __STDC__
 375 # include <stddef.h>
 376 # include <stdlib.h>
 377 # define HAS_STDLIB
 378 # ifdef vms
 379 #  define __ID__(a)a
 380 # endif
 381 #else
 382 # ifndef BSD
 383 #  ifndef __TURBOC__
 384 #   include <memory.h>
 385 #  endif
 386 # endif
 387 # ifdef hpux
 388 #  ifdef _INCLUDE__STDC__
 389 #   include <stddef.h>
 390 #   include <stdlib.h>
 391 #  endif
 392 # endif
 393 # include <sys/types.h>
 394 # if !defined(MSDOS) || defined(__TURBOC__)
 395 #  define __ID__(a)a
 396 # endif
 397 #endif
 398
 399 #ifdef __ID__
 400 # define __CAT__(a,b)__ID__(a)b
 401 #else
 402 # define __CAT__(a,b)a##b
 403 #endif
 404
 405
 406 #ifdef BSD
 407 # include <strings.h>
 408 # define memcpy(a,b,n) (bcopy(b,a,n),a)
 409 # define memcmp(a,b,n) bcmp(a,b,n)
 410 # define strchr(s,c) index(s,c)
 411 # define strrchr(s,c) rindex(s,c)
 412 #else
 413 # include <string.h>
 414 #endif
 415
 416 #include <ctype.h>
 417 #include <math.h>
 418 #include <setjmp.h>
 419 #include <assert.h>
 420
 421
 422 #ifdef vms
 423
 424 #define LACK_LABS
 425 #define LACK_MEMMOVE
 426 #define LACK_MEMCPY
 427
 428 #else
 429
 430 #define LACK_LABS       /* Undefine these if your library has these */
 431 #define LACK_MEMMOVE
 432
 433 #endif
 434
 435
 436 typedef struct __p2c_jmp_buf {
 437     struct __p2c_jmp_buf *next;
 438     jmp_buf jbuf;
 439 } __p2c_jmp_buf;
 440
 441
 442 /* Warning: The following will not work if setjmp is used simultaneously.
 443    This also violates the ANSI restriction about using vars after longjmp,
 444    but a typical implementation of longjmp will get it right anyway. */
 445
 446 #ifndef FAKE_TRY
 447 # define TRY(x)         do { __p2c_jmp_buf __try_jb;  \
 448                              __try_jb.next = __top_jb;  \
 449                              if (!setjmp((__top_jb = &__try_jb)->jbuf)) {
 450 # define RECOVER(x)     __top_jb = __try_jb.next; } else {
 451 # define RECOVER2(x,L)  __top_jb = __try_jb.next; } else {  \
 452                              if (0) { L: __top_jb = __try_jb.next; }
 453 # define ENDTRY(x)      } } while (0)
 454 #else
 455 # define TRY(x)         if (1) {
 456 # define RECOVER(x)     } else do {
 457 # define RECOVER2(x,L)  } else do { L: ;
 458 # define ENDTRY(x)      } while (0)
 459 #endif
 460
 461
 462
 463 #ifdef M_XENIX  /* avoid compiler bug */
 464 # define SHORT_MAX  (32767)
 465 # define SHORT_MIN  (-32768)
 466 #endif
 467
 468
 469 /* The following definitions work only on twos-complement machines */
 470 #ifndef SHORT_MAX
 471 # define SHORT_MAX  ((short)(((unsigned short) -1) >> 1))
 472 # define SHORT_MIN  (~SHORT_MAX)
 473 #endif
 474
 475 #ifndef INT_MAX
 476 # define INT_MAX    ((int)(((unsigned int) -1) >> 1))
 477 # define INT_MIN    (~INT_MAX)
 478 #endif
 479
 480 #ifndef LONG_MAX
 481 # define LONG_MAX   ((long)(((unsigned long) -1) >> 1))
 482 # define LONG_MIN   (~LONG_MAX)
 483 #endif
 484
 485 #ifndef SEEK_SET
 486 # define SEEK_SET   0
 487 # define SEEK_CUR   1
 488 # define SEEK_END   2
 489 #endif
 490
 491 #ifndef EXIT_SUCCESS
 492 # ifdef vms
 493 #  define EXIT_SUCCESS  1
 494 #  define EXIT_FAILURE  (02000000000L)
 495 # else
 496 #  define EXIT_SUCCESS  0
 497 #  define EXIT_FAILURE  1
 498 # endif
 499 #endif
 500
 501
 502 #define SETBITS  32
 503
 504
 505 #ifdef __STDC__
 506 # ifndef vms
 507 #  define Signed    signed
 508 # else
 509 #  define Signed
 510 # endif
 511 # define Void       void      /* Void f() = procedure */
 512 # ifndef Const
 513 #  define Const     const
 514 # endif
 515 # ifndef Volatile
 516 # define Volatile  volatile
 517 # endif
 518 # define PP(x)      x         /* function prototype */
 519 # define PV()       (void)    /* null function prototype */
 520 typedef void *Anyptr;
 521 #else
 522 # define Signed
 523 # define Void       void
 524 # ifndef Const
 525 #  define Const
 526 # endif
 527 # ifndef Volatile
 528 #  define Volatile
 529 # endif
 530 # define PP(x)      ()
 531 # define PV()       ()
 532 typedef char *Anyptr;
 533 #endif
 534
 535 #ifdef __GNUC__
 536 # define Inline     inline
 537 #else
 538 # define Inline
 539 #endif
 540
 541 #define Register    register  /* Register variables */
 542 #define Char        char      /* Characters (not bytes) */
 543
 544 #ifndef Static
 545 # define Static     static    /* Private global funcs and vars */
 546 #endif
 547
 548 #ifndef Local
 549 # define Local      static    /* Nested functions */
 550 #endif
 551
 552 typedef Signed   char schar;
 553 typedef unsigned char _uchar;
 554 typedef unsigned char boolean;
 555
 556 #ifndef true
 557 # define true    1
 558 # define false   0
 559 #endif
 560
 561
 562 typedef struct {
 563     Anyptr proc, link;
 564 } _PROCEDURE;
 565
 566 #ifndef _FNSIZE
 567 # define _FNSIZE  120
 568 #endif
 569
 570
 571 extern Void    PASCAL_MAIN  PP( (int, Char **) );
 572 extern Char    **P_argv;
 573 extern int     P_argc;
 574 extern short   P_escapecode;
 575 extern int     P_ioresult;
 576 extern __p2c_jmp_buf *__top_jb;
 577
 578
 579 #ifdef P2C_H_PROTO   /* if you have Ansi C but non-prototyped header files */
 580 extern Char    *strcat      PP( (Char *, Const Char *) );
 581 extern Char    *strchr      PP( (Const Char *, int) );
 582 extern int      strcmp      PP( (Const Char *, Const Char *) );
 583 extern Char    *strcpy      PP( (Char *, Const Char *) );
 584 extern size_t   strlen      PP( (Const Char *) );
 585 extern Char    *strncat     PP( (Char *, Const Char *, size_t) );
 586 extern int      strncmp     PP( (Const Char *, Const Char *, size_t) );
 587 extern Char    *strncpy     PP( (Char *, Const Char *, size_t) );
 588 extern Char    *strrchr     PP( (Const Char *, int) );
 589
 590 extern Anyptr   memchr      PP( (Const Anyptr, int, size_t) );
 591 extern Anyptr   memmove     PP( (Anyptr, Const Anyptr, size_t) );
 592 extern Anyptr   memset      PP( (Anyptr, int, size_t) );
 593 #ifndef memcpy
 594 extern Anyptr   memcpy      PP( (Anyptr, Const Anyptr, size_t) );
 595 extern int      memcmp      PP( (Const Anyptr, Const Anyptr, size_t) );
 596 #endif
 597
 598 extern int      atoi        PP( (Const Char *) );
 599 extern double   atof        PP( (Const Char *) );
 600 extern long     atol        PP( (Const Char *) );
 601 extern double   strtod      PP( (Const Char *, Char **) );
 602 extern long     strtol      PP( (Const Char *, Char **, int) );
 603 #endif /*P2C_H_PROTO*/
 604
 605 #define HAS_STDLIB
 606 #ifndef HAS_STDLIB
 607 extern Anyptr   malloc      PP( (size_t) );
 608 extern Void     free        PP( (Anyptr) );
 609 #endif
 610
 611 extern int      _OutMem     PV();
 612 extern int      _CaseCheck  PV();
 613 extern int      _NilCheck   PV();
 614 extern int      _Escape     PP( (int) );
 615 extern int      _EscIO      PP( (int) );
 616
 617 extern long     ipow        PP( (long, long) );
 618 extern Char    *strsub      PP( (Char *, Char *, int, int) );
 619 extern Char    *strltrim    PP( (Char *) );
 620 extern Char    *strrtrim    PP( (Char *) );
 621 extern Char    *strrpt      PP( (Char *, Char *, int) );
 622 extern Char    *strpad      PP( (Char *, Char *, int, int) );
 623 extern int      strpos2     PP( (Char *, Char *, int) );
 624 extern long     memavail    PV();
 625 extern int      P_peek      PP( (FILE *) );
 626 extern int      P_eof       PP( (FILE *) );
 627 extern int      P_eoln      PP( (FILE *) );
 628 extern Void     P_readpaoc  PP( (FILE *, Char *, int) );
 629 extern Void     P_readlnpaoc PP( (FILE *, Char *, int) );
 630 extern long     P_maxpos    PP( (FILE *) );
 631 extern Char    *P_trimname  PP( (Char *, int) );
 632 extern long    *P_setunion  PP( (long *, long *, long *) );
 633 extern long    *P_setint    PP( (long *, long *, long *) );
 634 extern long    *P_setdiff   PP( (long *, long *, long *) );
 635 extern long    *P_setxor    PP( (long *, long *, long *) );
 636 extern int      P_inset     PP( (unsigned, long *) );
 637 extern int      P_setequal  PP( (long *, long *) );
 638 extern int      P_subset    PP( (long *, long *) );
 639 extern long    *P_addset    PP( (long *, unsigned) );
 640 extern long    *P_addsetr   PP( (long *, unsigned, unsigned) );
 641 extern long    *P_remset    PP( (long *, unsigned) );
 642 extern long    *P_setcpy    PP( (long *, long *) );
 643 extern long    *P_expset    PP( (long *, long) );
 644 extern long     P_packset   PP( (long *) );
 645 extern int      P_getcmdline PP( (int l, int h, Char *line) );
 646 extern Void     TimeStamp   PP( (int *Day, int *Month, int *Year,
 647                                  int *Hour, int *Min, int *Sec) );
 648 extern Void     P_sun_argv  PP( (char *, int, int) );
 649
 650
 651 /* I/O error handling */
 652 #define _CHKIO(cond,ior,val,def)  ((cond) ? P_ioresult=0,(val)  \
 653                                           : P_ioresult=(ior),(def))
 654 #define _SETIO(cond,ior)          (P_ioresult = (cond) ? 0 : (ior))
 655
 656 /* Following defines are suitable for the HP Pascal operating system */
 657 #define FileNotFound     10
 658 #define FileNotOpen      13
 659 #define FileWriteError   38
 660 #define BadInputFormat   14
 661 #define EndOfFile        30
 662
 663 /* Creating temporary files */
 664 #if (defined(BSD) || defined(NO_TMPFILE)) && !defined(HAVE_TMPFILE)
 665 # define tmpfile()  (fopen(tmpnam(NULL), "w+"))
 666 #endif
 667
 668 /* File buffers */
 669 #define FILEBUF(f,sc,type) sc int __CAT__(f,_BFLAGS);   \
 670                            sc type __CAT__(f,_BUFFER)
 671 #define FILEBUFNC(f,type)  int __CAT__(f,_BFLAGS);   \
 672                            type __CAT__(f,_BUFFER)
 673
 674 #define RESETBUF(f,type)   (__CAT__(f,_BFLAGS) = 1)
 675 #define SETUPBUF(f,type)   (__CAT__(f,_BFLAGS) = 0)
 676
 677 #define GETFBUF(f,type)    (*((__CAT__(f,_BFLAGS) == 1 &&   \
 678                                ((__CAT__(f,_BFLAGS) = 2),   \
 679                                 fread(&__CAT__(f,_BUFFER),  \
 680                                       sizeof(type),1,(f)))),\
 681                               &__CAT__(f,_BUFFER)))
 682 #define AGETFBUF(f,type)   ((__CAT__(f,_BFLAGS) == 1 &&   \
 683                              ((__CAT__(f,_BFLAGS) = 2),   \
 684                               fread(__CAT__(f,_BUFFER),  \
 685                                     sizeof(type),1,(f)))),\
 686                             __CAT__(f,_BUFFER))
 687
 688 #define PUTFBUF(f,type,v)  (GETFBUF(f,type) = (v))
 689 #define CPUTFBUF(f,v)      (PUTFBUF(f,char,v))
 690 #define APUTFBUF(f,type,v) (memcpy(AGETFBUF(f,type), (v),  \
 691                                    sizeof(__CAT__(f,_BUFFER))))
 692
 693 #define GET(f,type)        (__CAT__(f,_BFLAGS) == 1 ?   \
 694                             fread(&__CAT__(f,_BUFFER),sizeof(type),1,(f)) :  \
 695                             (__CAT__(f,_BFLAGS) = 1))
 696
 697 #define PUT(f,type)        (fwrite(&__CAT__(f,_BUFFER),sizeof(type),1,(f)),  \
 698                             (__CAT__(f,_BFLAGS) = 0))
 699 #define CPUT(f)            (PUT(f,char))
 700
 701 #define BUFEOF(f)          (__CAT__(f,_BFLAGS) != 2 && P_eof(f))
 702 #define BUFFPOS(f)         (ftell(f) - (__CAT__(f,_BFLAGS) == 2))
 703
 704 typedef struct {
 705     FILE *f;
 706     FILEBUFNC(f,Char);
 707     Char name[_FNSIZE];
 708 } _TEXT;
 709
 710 /* Memory allocation */
 711 #ifdef __GCC__
 712 # define Malloc(n)  (malloc(n) ?: (Anyptr)_OutMem())
 713 #else
 714 extern Anyptr __MallocTemp__;
 715 # define Malloc(n)  ((__MallocTemp__ = malloc(n)) ? __MallocTemp__ : (Anyptr)_OutMem())
 716 #endif
 717 #define FreeR(p)    (free((Anyptr)(p)))    /* used if arg is an rvalue */
 718 #define Free(p)     (free((Anyptr)(p)), (p)=NULL)
 719
 720 /* sign extension */
 721 #define SEXT(x,n)   ((x) | -(((x) & (1L<<((n)-1))) << 1))
 722
 723 /* packed arrays */   /* BEWARE: these are untested! */
 724 #define P_getbits_UB(a,i,n,L)   ((int)((a)[(i)>>(L)-(n)] >>   \
 725                                        (((~(i))&((1<<(L)-(n))-1)) << (n)) &  \
 726                                        (1<<(1<<(n)))-1))
 727
 728 #define P_getbits_SB(a,i,n,L)   ((int)((a)[(i)>>(L)-(n)] <<   \
 729                                        (16 - ((((~(i))&((1<<(L)-(n))-1))+1) <<\
 730                                               (n)) >> (16-(1<<(n))))))
 731
 732 #define P_putbits_UB(a,i,x,n,L) ((a)[(i)>>(L)-(n)] |=   \
 733                                  (x) << (((~(i))&((1<<(L)-(n))-1)) << (n)))
 734
 735 #define P_putbits_SB(a,i,x,n,L) ((a)[(i)>>(L)-(n)] |=   \
 736                                  ((x) & (1<<(1<<(n)))-1) <<   \
 737                                  (((~(i))&((1<<(L)-(n))-1)) << (n)))
 738
 739 #define P_clrbits_B(a,i,n,L)    ((a)[(i)>>(L)-(n)] &=   \
 740                                  ~( ((1<<(1<<(n)))-1) <<   \
 741                                    (((~(i))&((1<<(L)-(n))-1)) << (n))) )
 742
 743 /* small packed arrays */
 744 #define P_getbits_US(v,i,n)     ((int)((v) >> ((i)<<(n)) & (1<<(1<<(n)))-1))
 745 #define P_getbits_SS(v,i,n)     ((int)((long)(v) << (SETBITS - (((i)+1) << (n))) >> (SETBITS-(1<<(n)))))
 746 #define P_putbits_US(v,i,x,n)   ((v) |= (x) << ((i) << (n)))
 747 #define P_putbits_SS(v,i,x,n)   ((v) |= ((x) & (1<<(1<<(n)))-1) << ((i)<<(n)))
 748 #define P_clrbits_S(v,i,n)      ((v) &= ~( ((1<<(1<<(n)))-1) << ((i)<<(n)) ))
 749
 750 #define P_max(a,b)   ((a) > (b) ? (a) : (b))
 751 #define P_min(a,b)   ((a) < (b) ? (a) : (b))
 752
 753
 754 /* Fix ANSI-isms */
 755
 756 #ifdef LACK_LABS
 757 # ifndef labs
 758 #  define labs  my_labs
 759    extern long my_labs PP( (long) );
 760 # endif
 761 #endif
 762
 763 #ifdef LACK_MEMMOVE
 764 # ifndef memmove
 765 #  define memmove  my_memmove
 766    extern Anyptr my_memmove PP( (Anyptr, Const Anyptr, size_t) );
 767 # endif
 768 #endif
 769
 770 #ifdef LACK_MEMCPY
 771 # ifndef memcpy
 772 #  define memcpy  my_memcpy
 773    extern Anyptr my_memcpy PP( (Anyptr, Const Anyptr, size_t) );
 774 # endif
 775 # ifndef memcmp
 776 #  define memcmp  my_memcmp
 777    extern int my_memcmp PP( (Const Anyptr, Const Anyptr, size_t) );
 778 # endif
 779 # ifndef memset
 780 #  define memset  my_memset
 781    extern Anyptr my_memset PP( (Anyptr, int, size_t) );
 782 # endif
 783 #endif
 784
 785 /* Fix toupper/tolower on Suns and other stupid BSD systems */
 786 #ifdef toupper
 787 # undef toupper
 788 # undef tolower
 789 # define toupper(c)   my_toupper(c)
 790 # define tolower(c)   my_tolower(c)
 791 #endif
 792
 793 #ifndef _toupper
 794 # if 'A' == 65 && 'a' == 97
 795 #  define _toupper(c)  ((c)-'a'+'A')
 796 #  define _tolower(c)  ((c)-'A'+'a')
 797 # else
 798 #  ifdef toupper
 799 #   undef toupper   /* hope these are shadowing real functions, */
 800 #   undef tolower   /* because my_toupper calls _toupper! */
 801 #  endif
 802 #  define _toupper(c)  toupper(c)
 803 #  define _tolower(c)  tolower(c)
 804 # endif
 805 #endif
 806
 807
 808 #endif    /* P2C_H */
 809
 810
 811
 812 /* End. */
 813
 814
 815
 816
 817 /************************************************************
 818 *                                                           *
 819 *          p2clib.c                                           *
 820 *                                                           *
 821 ************************************************************/
 822
 823 /* Run-time library for use with "p2c", the Pascal to C translator */
 824
 825 /* "p2c"  Copyright (C) 1989, 1990, 1991 Free Software Foundation.
 826  * By Dave Gillespie, daveg@csvax.cs.caltech.edu.  Version --VERSION--.
 827  * This file may be copied, modified, etc. in any way.  It is not restricted
 828  * by the licence agreement accompanying p2c itself.
 829  */
 830
 831
 832
 833  /* #include "p2c.h" */
 834
 835
 836 #ifndef NO_TIME
 837 # include <time.h>
 838 #endif
 839
 840
 841 #define Isspace(c)  isspace(c)      /* or "((c) == ' ')" if preferred */
 842
 843
 844
 845
 846 int P_argc;
 847 char **P_argv;
 848
 849 short P_escapecode;
 850 int P_ioresult;
 851
 852 long EXCP_LINE;    /* Used by Pascal workstation system */
 853
 854 Anyptr __MallocTemp__;
 855
 856 __p2c_jmp_buf *__top_jb;
 857
 858
 859
 860
 861 void PASCAL_MAIN(argc, argv)
 862 int argc;
 863 char **argv;
 864 {
 865     P_argc = argc;
 866     P_argv = argv;
 867     __top_jb = NULL;
 868
 869 #ifdef LOCAL_INIT
 870     LOCAL_INIT();
 871 #endif
 872 }
 873
 874
 875
 876
 877
 878 /* In case your system lacks these... */
 879
 880 long my_labs(x)
 881 long x;
 882 {
 883     return((x > 0) ? x : -x);
 884 }
 885
 886
 887 #ifdef __STDC__
 888 Anyptr my_memmove(Anyptr d, Const Anyptr s, size_t n)
 889 #else
 890 Anyptr my_memmove(d, s, n)
 891 Anyptr d, s;
 892 register long n;
 893 #endif
 894 {
 895     register char *dd = (char *)d, *ss = (char *)s;
 896     if (dd < ss || dd - ss >= n) {
 897         memcpy(dd, ss, n);
 898     } else if (n > 0) {
 899         dd += n;
 900         ss += n;
 901         while (--n >= 0)
 902             *--dd = *--ss;
 903     }
 904     return d;
 905 }
 906
 907
 908 #ifdef __STDC__
 909 Anyptr my_memcpy(Anyptr d, Const Anyptr s, size_t n)
 910 #else
 911 Anyptr my_memcpy(d, s, n)
 912 Anyptr d, s;
 913 register long n;
 914 #endif
 915 {
 916     register char *ss = (char *)s, *dd = (char *)d;
 917     while (--n >= 0)
 918         *dd++ = *ss++;
 919     return d;
 920 }
 921
 922 #ifdef __STDC__
 923 int my_memcmp(Const Anyptr s1, Const Anyptr s2, size_t n)
 924 #else
 925 int my_memcmp(s1, s2, n)
 926 Anyptr s1, s2;
 927 register long n;
 928 #endif
 929 {
 930     register char *a = (char *)s1, *b = (char *)s2;
 931     register int i;
 932     while (--n >= 0)
 933         if ((i = (*a++) - (*b++)) != 0)
 934             return i;
 935     return 0;
 936 }
 937
 938 #ifdef __STDC__
 939 Anyptr my_memset(Anyptr d, int c, size_t n)
 940 #else
 941 Anyptr my_memset(d, c, n)
 942 Anyptr d;
 943 register int c;
 944 register long n;
 945 #endif
 946 {
 947     register char *dd = (char *)d;
 948     while (--n >= 0)
 949         *dd++ = c;
 950     return d;
 951 }
 952
 953
 954 int my_toupper(c)
 955 int c;
 956 {
 957     if (islower(c))
 958         return _toupper(c);
 959     else
 960         return c;
 961 }
 962
 963
 964 int my_tolower(c)
 965 int c;
 966 {
 967     if (isupper(c))
 968         return _tolower(c);
 969     else
 970         return c;
 971 }
 972
 973
 974
 975
 976 long ipow(a, b)
 977 long a, b;
 978 {
 979     long v;
 980
 981     if (a == 0 || a == 1)
 982         return a;
 983     if (a == -1)
 984         return (b & 1) ? -1 : 1;
 985     if (b < 0)
 986         return 0;
 987     if (a == 2)
 988         return 1 << b;
 989     v = (b & 1) ? a : 1;
 990     while ((b >>= 1) > 0) {
 991         a *= a;
 992         if (b & 1)
 993             v *= a;
 994     }
 995     return v;
 996 }
 997
 998
 999
1000
1001 /* Common string functions: */
1002
1003 /* Store in "ret" the substring of length "len" starting from "pos" (1-based).
1004    Store a shorter or null string if out-of-range.  Return "ret". */
1005
1006 char *strsub(ret, s, pos, len)
1007 register char *ret, *s;
1008 register int pos, len;
1009 {
1010     register char *s2;
1011
1012     if (--pos < 0 || len <= 0) {
1013         *ret = 0;
1014         return ret;
1015     }
1016     while (pos > 0) {
1017         if (!*s++) {
1018             *ret = 0;
1019             return ret;
1020         }
1021         pos--;
1022     }
1023     s2 = ret;
1024     while (--len >= 0) {
1025         if (!(*s2++ = *s++))
1026             return ret;
1027     }
1028     *s2 = 0;
1029     return ret;
1030 }
1031
1032
1033 /* Return the index of the first occurrence of "pat" as a substring of "s",
1034    starting at index "pos" (1-based).  Result is 1-based, 0 if not found. */
1035
1036 int strpos2(s, pat, pos)
1037 char *s;
1038 register char *pat;
1039 register int pos;
1040 {
1041     register char *cp, ch;
1042     register int slen;
1043
1044     if (--pos < 0)
1045         return 0;
1046     slen = strlen(s) - pos;
1047     cp = s + pos;
1048     if (!(ch = *pat++))
1049         return 0;
1050     pos = strlen(pat);
1051     slen -= pos;
1052     while (--slen >= 0) {
1053         if (*cp++ == ch && !strncmp(cp, pat, pos))
1054             return cp - s;
1055     }
1056     return 0;
1057 }
1058
1059
1060 /* Case-insensitive version of strcmp. */
1061
1062 int strcicmp(s1, s2)
1063 register char *s1, *s2;
1064 {
1065     register unsigned char c1, c2;
1066
1067     while (*s1) {
1068         if (*s1++ != *s2++) {
1069             if (!s2[-1])
1070                 return 1;
1071             c1 = toupper(s1[-1]);
1072             c2 = toupper(s2[-1]);
1073             if (c1 != c2)
1074                 return c1 - c2;
1075         }
1076     }
1077     if (*s2)
1078         return -1;
1079     return 0;
1080 }
1081
1082
1083
1084
1085 /* HP and Turbo Pascal string functions: */
1086
1087 /* Trim blanks at left end of string. */
1088
1089 char *strltrim(s)
1090 register char *s;
1091 {
1092     while (Isspace(*s++)) ;
1093     return s - 1;
1094 }
1095
1096
1097 /* Trim blanks at right end of string. */
1098
1099 char *strrtrim(s)
1100 register char *s;
1101 {
1102     register char *s2 = s;
1103
1104     if (!*s)
1105         return s;
1106     while (*++s2) ;
1107     while (s2 > s && Isspace(*--s2))
1108         *s2 = 0;
1109     return s;
1110 }
1111
1112
1113 /* Store in "ret" "num" copies of string "s".  Return "ret". */
1114
1115 char *strrpt(ret, s, num)
1116 char *ret;
1117 register char *s;
1118 register int num;
1119 {
1120     register char *s2 = ret;
1121     register char *s1;
1122
1123     while (--num >= 0) {
1124         s1 = s;
1125         while ((*s2++ = *s1++)) ;
1126         s2--;
1127     }
1128     return ret;
1129 }
1130
1131
1132 /* Store in "ret" string "s" with enough pad chars added to reach "size". */
1133
1134 char *strpad(ret, s, padchar, num)
1135 char *ret;
1136 register char *s;
1137 register int padchar, num;
1138 {
1139     register char *d = ret;
1140
1141     if (s == d) {
1142         while (*d++) ;
1143     } else {
1144         while ((*d++ = *s++)) ;
1145     }
1146     num -= (--d - ret);
1147     while (--num >= 0)
1148         *d++ = padchar;
1149     *d = 0;
1150     return ret;
1151 }
1152
1153
1154 /* Copy the substring of length "len" from index "spos" of "s" (1-based)
1155    to index "dpos" of "d", lengthening "d" if necessary.  Length and
1156    indices must be in-range. */
1157
1158 void strmove(len, s, spos, d, dpos)
1159 register char *s, *d;
1160 register int len, spos, dpos;
1161 {
1162     s += spos - 1;
1163     d += dpos - 1;
1164     while (*d && --len >= 0)
1165         *d++ = *s++;
1166     if (len > 0) {
1167         while (--len >= 0)
1168             *d++ = *s++;
1169         *d = 0;
1170     }
1171 }
1172
1173
1174 /* Delete the substring of length "len" at index "pos" from "s".
1175    Delete less if out-of-range. */
1176
1177 void strdelete(s, pos, len)
1178 register char *s;
1179 register int pos, len;
1180 {
1181     register int slen;
1182
1183     if (--pos < 0)
1184         return;
1185     slen = strlen(s) - pos;
1186     if (slen <= 0)
1187         return;
1188     s += pos;
1189     if (slen <= len) {
1190         *s = 0;
1191         return;
1192     }
1193     while ((*s = s[len])) s++;
1194 }
1195
1196
1197 /* Insert string "src" at index "pos" of "dst". */
1198
1199 void strinsert(src, dst, pos)
1200 register char *src, *dst;
1201 register int pos;
1202 {
1203     register int slen, dlen;
1204
1205     if (--pos < 0)
1206         return;
1207     dlen = strlen(dst);
1208     dst += dlen;
1209     dlen -= pos;
1210     if (dlen <= 0) {
1211         strcpy(dst, src);
1212         return;
1213     }
1214     slen = strlen(src);
1215     do {
1216         dst[slen] = *dst;
1217         --dst;
1218     } while (--dlen >= 0);
1219     dst++;
1220     while (--slen >= 0)
1221         *dst++ = *src++;
1222 }
1223
1224
1225
1226
1227 /* File functions */
1228
1229 /* Peek at next character of input stream; return EOF at end-of-file. */
1230
1231 int P_peek(f)
1232 FILE *f;
1233 {
1234     int ch;
1235
1236     ch = getc(f);
1237     if (ch == EOF)
1238         return EOF;
1239     ungetc(ch, f);
1240     return (ch == '\n') ? ' ' : ch;
1241 }
1242
1243
1244 /* Check if at end of file, using Pascal "eof" semantics.  End-of-file for
1245    stdin is broken; remove the special case for it to be broken in a
1246    different way. */
1247
1248 int P_eof(f)
1249 FILE *f;
1250 {
1251     register int ch;
1252
1253     if (feof(f))
1254         return 1;
1255     if (f == stdin)
1256         return 0;    /* not safe to look-ahead on the keyboard! */
1257     ch = getc(f);
1258     if (ch == EOF)
1259         return 1;
1260     ungetc(ch, f);
1261     return 0;
1262 }
1263
1264
1265 /* Check if at end of line (or end of entire file). */
1266
1267 int P_eoln(f)
1268 FILE *f;
1269 {
1270     register int ch;
1271
1272     ch = getc(f);
1273     if (ch == EOF)
1274         return 1;
1275     ungetc(ch, f);
1276     return (ch == '\n');
1277 }
1278
1279
1280 /* Read a packed array of characters from a file. */
1281
1282 Void P_readpaoc(f, s, len)
1283 FILE *f;
1284 char *s;
1285 int len;
1286 {
1287     int ch;
1288
1289     for (;;) {
1290         if (len <= 0)
1291             return;
1292         ch = getc(f);
1293         if (ch == EOF || ch == '\n')
1294             break;
1295         *s++ = ch;
1296         --len;
1297     }
1298     while (--len >= 0)
1299         *s++ = ' ';
1300     if (ch != EOF)
1301         ungetc(ch, f);
1302 }
1303
1304 Void P_readlnpaoc(f, s, len)
1305 FILE *f;
1306 char *s;
1307 int len;
1308 {
1309     int ch;
1310
1311     for (;;) {
1312         ch = getc(f);
1313         if (ch == EOF || ch == '\n')
1314             break;
1315         if (len > 0) {
1316             *s++ = ch;
1317             --len;
1318         }
1319     }
1320     while (--len >= 0)
1321         *s++ = ' ';
1322 }
1323
1324
1325 /* Compute maximum legal "seek" index in file (0-based). */
1326
1327 long P_maxpos(f)
1328 FILE *f;
1329 {
1330     long savepos = ftell(f);
1331     long val;
1332
1333     if (fseek(f, 0L, SEEK_END))
1334         return -1;
1335     val = ftell(f);
1336     if (fseek(f, savepos, SEEK_SET))
1337         return -1;
1338     return val;
1339 }
1340
1341
1342 /* Use packed array of char for a file name. */
1343
1344 Char *P_trimname(fn, len)
1345 register Char *fn;
1346 register int len;
1347 {
1348     static Char fnbuf[256];
1349     register Char *cp = fnbuf;
1350
1351     while (--len >= 0 && *fn && !isspace(*fn))
1352         *cp++ = *fn++;
1353     *cp = 0;
1354     return fnbuf;
1355 }
1356
1357
1358
1359
1360
1361
1362 /* Sets are stored as an array of longs.  S[0] is the size of the set;
1363    S[N] is the N'th 32-bit chunk of the set.  S[0] equals the maximum
1364    I such that S[I] is nonzero.  S[0] is zero for an empty set.  Within
1365    each long, bits are packed from lsb to msb.  The first bit of the
1366    set is the element with ordinal value 0.  (Thus, for a "set of 5..99",
1367    the lowest five bits of the first long are unused and always zero.) */
1368
1369 /* (Sets with 32 or fewer elements are normally stored as plain longs.) */
1370
1371 long *P_setunion(d, s1, s2)         /* d := s1 + s2 */
1372 register long *d, *s1, *s2;
1373 {
1374     long *dbase = d++;
1375     register int sz1 = *s1++, sz2 = *s2++;
1376     while (sz1 > 0 && sz2 > 0) {
1377         *d++ = *s1++ | *s2++;
1378         sz1--, sz2--;
1379     }
1380     while (--sz1 >= 0)
1381         *d++ = *s1++;
1382     while (--sz2 >= 0)
1383         *d++ = *s2++;
1384     *dbase = d - dbase - 1;
1385     return dbase;
1386 }
1387
1388
1389 long *P_setint(d, s1, s2)           /* d := s1 * s2 */
1390 register long *d, *s1, *s2;
1391 {
1392     long *dbase = d++;
1393     register int sz1 = *s1++, sz2 = *s2++;
1394     while (--sz1 >= 0 && --sz2 >= 0)
1395         *d++ = *s1++ & *s2++;
1396     while (--d > dbase && !*d) ;
1397     *dbase = d - dbase;
1398     return dbase;
1399 }
1400
1401
1402 long *P_setdiff(d, s1, s2)          /* d := s1 - s2 */
1403 register long *d, *s1, *s2;
1404 {
1405     long *dbase = d++;
1406     register int sz1 = *s1++, sz2 = *s2++;
1407     while (--sz1 >= 0 && --sz2 >= 0)
1408         *d++ = *s1++ & ~*s2++;
1409     if (sz1 >= 0) {
1410         while (sz1-- >= 0)
1411             *d++ = *s1++;
1412     }
1413     while (--d > dbase && !*d) ;
1414     *dbase = d - dbase;
1415     return dbase;
1416 }
1417
1418
1419 long *P_setxor(d, s1, s2)         /* d := s1 / s2 */
1420 register long *d, *s1, *s2;
1421 {
1422     long *dbase = d++;
1423     register int sz1 = *s1++, sz2 = *s2++;
1424     while (sz1 > 0 && sz2 > 0) {
1425         *d++ = *s1++ ^ *s2++;
1426         sz1--, sz2--;
1427     }
1428     while (--sz1 >= 0)
1429         *d++ = *s1++;
1430     while (--sz2 >= 0)
1431         *d++ = *s2++;
1432     while (--d > dbase && !*d) ;
1433     *dbase = d - dbase;
1434     return dbase;
1435 }
1436
1437
1438 int P_inset(val, s)                 /* val IN s */
1439 register unsigned val;
1440 register long *s;
1441 {
1442     register int bit;
1443     bit = val % SETBITS;
1444     val /= SETBITS;
1445     if (val < *s++ && ((1<<bit) & s[val]))
1446         return 1;
1447     return 0;
1448 }
1449
1450
1451 long *P_addset(s, val)              /* s := s + [val] */
1452 register long *s;
1453 register unsigned val;
1454 {
1455     register long *sbase = s;
1456     register int bit, size;
1457     bit = val % SETBITS;
1458     val /= SETBITS;
1459     size = *s;
1460     if (++val > size) {
1461         s += size;
1462         while (val > size)
1463             *++s = 0, size++;
1464         *sbase = size;
1465     } else
1466         s += val;
1467     *s |= 1<<bit;
1468     return sbase;
1469 }
1470
1471
1472 long *P_addsetr(s, v1, v2)              /* s := s + [v1..v2] */
1473 register long *s;
1474 register unsigned v1, v2;
1475 {
1476     register long *sbase = s;
1477     register int b1, b2, size;
1478     if ((int)v1 > (int)v2)
1479         return sbase;
1480     b1 = v1 % SETBITS;
1481     v1 /= SETBITS;
1482     b2 = v2 % SETBITS;
1483     v2 /= SETBITS;
1484     size = *s;
1485     v1++;
1486     if (++v2 > size) {
1487         while (v2 > size)
1488             s[++size] = 0;
1489         s[v2] = 0;
1490         *s = v2;
1491     }
1492     s += v1;
1493     if (v1 == v2) {
1494         *s |= (~((-2)<<(b2-b1))) << b1;
1495     } else {
1496         *s++ |= (-1) << b1;
1497         while (++v1 < v2)
1498             *s++ = -1;
1499         *s |= ~((-2) << b2);
1500     }
1501     return sbase;
1502 }
1503
1504
1505 long *P_remset(s, val)              /* s := s - [val] */
1506 register long *s;
1507 register unsigned val;
1508 {
1509     register int bit;
1510     bit = val % SETBITS;
1511     val /= SETBITS;
1512     if (++val <= *s) {
1513         if (!(s[val] &= ~(1<<bit)))
1514             while (*s && !s[*s])
1515                 (*s)--;
1516     }
1517     return s;
1518 }
1519
1520
1521 int P_setequal(s1, s2)              /* s1 = s2 */
1522 register long *s1, *s2;
1523 {
1524     register int size = *s1++;
1525     if (*s2++ != size)
1526         return 0;
1527     while (--size >= 0) {
1528         if (*s1++ != *s2++)
1529             return 0;
1530     }
1531     return 1;
1532 }
1533
1534
1535 int P_subset(s1, s2)                /* s1 <= s2 */
1536 register long *s1, *s2;
1537 {
1538     register int sz1 = *s1++, sz2 = *s2++;
1539     if (sz1 > sz2)
1540         return 0;
1541     while (--sz1 >= 0) {
1542         if (*s1++ & ~*s2++)
1543             return 0;
1544     }
1545     return 1;
1546 }
1547
1548
1549 long *P_setcpy(d, s)                /* d := s */
1550 register long *d, *s;
1551 {
1552     register long *save_d = d;
1553
1554 #ifdef SETCPY_MEMCPY
1555     memcpy(d, s, (*s + 1) * sizeof(long));
1556 #else
1557     register int i = *s + 1;
1558     while (--i >= 0)
1559         *d++ = *s++;
1560 #endif
1561     return save_d;
1562 }
1563
1564
1565 /* s is a "smallset", i.e., a 32-bit or less set stored
1566    directly in a long. */
1567
1568 long *P_expset(d, s)                /* d := s */
1569 register long *d;
1570 register long s;
1571 {
1572     if (s) {
1573         d[1] = s;
1574         *d = 1;
1575     } else
1576         *d = 0;
1577     return d;
1578 }
1579
1580
1581 long P_packset(s)                   /* convert s to a small-set */
1582 register long *s;
1583 {
1584     if (*s++)
1585         return *s;
1586     else
1587         return 0;
1588 }
1589
1590
1591
1592
1593
1594 /* Oregon Software Pascal extensions, courtesy of William Bader */
1595
1596 int P_getcmdline(l, h, line)
1597 int l, h;
1598 Char *line;
1599 {
1600     int i, len;
1601     char *s;
1602
1603     h = h - l + 1;
1604     len = 0;
1605     for(i = 1; i < P_argc; i++) {
1606         s = P_argv[i];
1607         while (*s) {
1608             if (len >= h) return len;
1609             line[len++] = *s++;
1610         }
1611         if (len >= h) return len;
1612         line[len++] = ' ';
1613     }
1614     return len;
1615 }
1616
1617 Void TimeStamp(Day, Month, Year, Hour, Min, Sec)
1618 int *Day, *Month, *Year, *Hour, *Min, *Sec;
1619 {
1620 #ifndef NO_TIME
1621     struct tm *tm;
1622     long clock;
1623
1624     time(&clock);
1625     tm = localtime(&clock);
1626     *Day = tm->tm_mday;
1627     *Month = tm->tm_mon + 1;            /* Jan = 0 */
1628     *Year = tm->tm_year;
1629     if (*Year < 1900)
1630         *Year += 1900;     /* year since 1900 */
1631     *Hour = tm->tm_hour;
1632     *Min = tm->tm_min;
1633     *Sec = tm->tm_sec;
1634 #endif
1635 }
1636
1637
1638
1639
1640 /* SUN Berkeley Pascal extensions */
1641
1642 Void P_sun_argv(s, len, n)
1643 register char *s;
1644 register int len, n;
1645 {
1646     register char *cp;
1647
1648     if ((unsigned)n < P_argc)
1649       cp = P_argv[n];
1650     else
1651       cp = "";
1652     while (*cp && --len >= 0)
1653       *s++ = *cp++;
1654     while (--len >= 0)
1655       *s++ = ' ';
1656 }
1657
1658
1659
1660
1661 int _OutMem()
1662 {
1663     return _Escape(-2);
1664 }
1665
1666 int _CaseCheck()
1667 {
1668     return _Escape(-9);
1669 }
1670
1671 int _NilCheck()
1672 {
1673     return _Escape(-3);
1674 }
1675
1676
1677
1678
1679
1680 /* The following is suitable for the HP Pascal operating system.
1681    It might want to be revised when emulating another system. */
1682
1683 char *_ShowEscape(buf, code, ior, prefix)
1684 char *buf, *prefix;
1685 int code, ior;
1686 {
1687     char *bufp;
1688
1689     if (prefix && *prefix) {
1690         strcpy(buf, prefix);
1691         strcat(buf, ": ");
1692         bufp = buf + strlen(buf);
1693     } else {
1694         bufp = buf;
1695     }
1696     if (code == -10) {
1697         sprintf(bufp, "Pascal system I/O error %d", ior);
1698         switch (ior) {
1699             case 3:
1700                 strcat(buf, " (illegal I/O request)");
1701                 break;
1702             case 7:
1703                 strcat(buf, " (bad file name)");
1704                 break;
1705             case FileNotFound:   /*10*/
1706                 strcat(buf, " (file not found)");
1707                 break;
1708             case FileNotOpen:    /*13*/
1709                 strcat(buf, " (file not open)");
1710                 break;
1711             case BadInputFormat: /*14*/
1712                 strcat(buf, " (bad input format)");
1713                 break;
1714             case 24:
1715                 strcat(buf, " (not open for reading)");
1716                 break;
1717             case 25:
1718                 strcat(buf, " (not open for writing)");
1719                 break;
1720             case 26:
1721                 strcat(buf, " (not open for direct access)");
1722                 break;
1723             case 28:
1724                 strcat(buf, " (string subscript out of range)");
1725                 break;
1726             case EndOfFile:      /*30*/
1727                 strcat(buf, " (end-of-file)");
1728                 break;
1729             case FileWriteError: /*38*/
1730                 strcat(buf, " (file write error)");
1731                 break;
1732         }
1733     } else {
1734         sprintf(bufp, "Pascal system error %d", code);
1735         switch (code) {
1736             case -2:
1737                 strcat(buf, " (out of memory)");
1738                 break;
1739             case -3:
1740                 strcat(buf, " (reference to NIL pointer)");
1741                 break;
1742             case -4:
1743                 strcat(buf, " (integer overflow)");
1744                 break;
1745             case -5:
1746                 strcat(buf, " (divide by zero)");
1747                 break;
1748             case -6:
1749                 strcat(buf, " (real math overflow)");
1750                 break;
1751             case -8:
1752                 strcat(buf, " (value range error)");
1753                 break;
1754             case -9:
1755                 strcat(buf, " (CASE value range error)");
1756                 break;
1757             case -12:
1758                 strcat(buf, " (bus error)");
1759                 break;
1760             case -20:
1761                 strcat(buf, " (stopped by user)");
1762                 break;
1763         }
1764     }
1765     return buf;
1766 }
1767
1768
1769 int _Escape(code)
1770 int code;
1771 {
1772     char buf[100];
1773
1774     P_escapecode = code;
1775     if (__top_jb) {
1776         __p2c_jmp_buf *jb = __top_jb;
1777         __top_jb = jb->next;
1778         longjmp(jb->jbuf, 1);
1779     }
1780     if (code == 0)
1781         exit(EXIT_SUCCESS);
1782     if (code == -1)
1783         exit(EXIT_FAILURE);
1784     fprintf(stderr, "%s\n", _ShowEscape(buf, P_escapecode, P_ioresult, ""));
1785     exit(EXIT_FAILURE);
1786 }
1787
1788 int _EscIO(code)
1789 int code;
1790 {
1791     P_ioresult = code;
1792     return _Escape(-10);
1793 }
1794
1795
1796
1797
1798 /* End. */
1799
1800
1801
1802
1803
1804 /************************************************************
1805 *                                                           *
1806 *          date.c                                           *
1807 *                                                           *
1808 ************************************************************/
1809
1810
1811
1812 /*#include <time.h>
1813 #include <stdio.h>*/
1814 static  char *mon[]={"JAN","FEB","MAR","APR","MAY","JUN",
1815                 "JUL","AUG","SEP","OCT","NOV","DEC"};
1816 /* PROCEDURE DATE(VAR DATESTRING:PACKED ARRAY[1..11] OF CHAR);EXTERN; */
1817 /* activate DATE by removing comment brackets if necessary */
1818 /***/
1819
1820 void Date(string)
1821 char* string;
1822 {
1823   time_t tt;
1824   struct tm *t;
1825   time(&tt);
1826   t=localtime(&tt);
1827   sprintf(string,"%d-%s-19%d\n",t->tm_mday,mon[t->tm_mon],t->tm_year);
1828 }
1829
1830
1831 /* p2c: dssp.p, line 295:
1832  * Note: Unexpected name "tapein" in program header [262] */
1833 /* p2c: dssp.p, line 295:
1834  * Note: Unexpected name "tapeout" in program header [262] */
1835
1836
1837 /*--------------------------------------------------------------------*/
1838 /* PROGRAM FATAL ERROR EXIT LABEL */
1839 /*******************  MATHEMATICAL CONSTANTS  **************************
1840  YVERTEX, - ARE Y,Z-COMPONENTS OF THE FIRST ICOSAHEDRON VERTEX. THE
1841  ZVERTEX    X-COMPONENT IS 0.
1842  EPS      - NUMERICAL TOLERANCE
1843   --------------------------------------------------------------------*/
1844
1845 #define PIHALF          1.570796
1846 #define PI              3.141593
1847 #define TWOPI           6.283185
1848 #define FOURPI          12.56637
1849 #define RADIAN          57.29578
1850 #define YVERTEX         0.8506508
1851 #define ZVERTEX         0.5257311
1852 #define EPS             0.00001
1853 /***/
1854 /***************  ARRAY DIMENSIONING CONSTANTS  ***********************
1855  NMAX     - MAXIMUM NUMBER OF AMINOACID RESIDUES IN ARRAY CHAIN
1856  MAXATOM  - MAXIMUM NUMBER OF SIDECHAIN ATOMS IN ARRAY SIDECHAIN
1857  MAXBRIDGE- MAXIMUM NUMBER OF BRIDGES IN ARRAY BRIDGETABLE
1858  NFACE,   - NUMBER OF FACES OF POLYHEDRON. THE COORDINATES OF THE CENTRE
1859  ORDER      OF EACH TRIANGULAR FACE ARE STORED IN ARRAY P, THE AREA
1860              IS STORED IN ARRAY WP IN PROCEDURE FLAGACCESS. NFACE MUST
1861              BE OF THE FORM NFACE=20*(4**ORDER), ORDER=0,1,2,...
1862              THE ACCURACY OF THE SOLVENT ACCESSIBLE SURFACE OF EACH
1863              AMINOACID RESIDUE IS ONE ANGSTROM**2 FOR ORDER=2,NFACE=320.
1864  MAXPACK  - MAXIMUM NUMBER OF PROTEIN ATOMS WHICH CAN INTRUDE INTO
1865              SOLVENT AROUND ANY GIVEN TEST ATOM. THE COORDINATES OF
1866              THESE ATOMS ARE STORED IN ARRAY X, THEIR RADII IN ARRAY RX
1867              IN PROCEDURE SURFACE.
1868  MAXHIST  - NUMBER OF SLOTS IN ARRAYS HELIXHIST AND BETAHIST USED FOR
1869              LENGTH STATISTICS OF SECONDARY STRUCTURE.
1870  MAXSS    - MAXIMUM NUMBER OF SSBOND RECORDS ON INPUT FILE. THE
1871              DISULFIDE BOND ARE SAVED IN ARRAY SSBONDS.
1872   --------------------------------------------------------------------*/
1873
1874 #define NMAX            6000
1875
1876 #define MAXATOM         40000L
1877
1878 #define MAXBRIDGE       300
1879 #define NFACE           320
1880 #define ORDER           2
1881 #define MAXPACK         200
1882 #define MAXHIST         30
1883 #define MAXSS           100
1884 /***/
1885 /*********************  PHYSICAL CONSTANTS   **************************
1886  RN       - RADIUS OF PEPTIDE NITROGEN ATOM
1887  RCA      - RADIUS OF PEPTIDE ALPHA-CARBON ATOM
1888  RC       - RADIUS OF PEPTIDE C'-CARBON ATOM
1889  RO       - RADIUS OF PEPTIDE OXYGEN ATOM
1890  RSIDEATOM- RADIUS OF SIDECHAIN ATOM
1891  RWATER   - RADIUS OF WATER MOLECULE
1892  SSDIST   - MAXIMUM ALLOWED DISTANCE OF DISULFIDE BRIDGE
1893  BREAKDIST- MAXIMUM ALLOWED PEPTIDE BOND LENGTH. IF DISTANCE IS
1894              GREATER A POLYPEPTIDE CHAIN INTERRUPTION IS ASSUMED.
1895  RESRAD   - MAXIMUM RADIUS OF A SPHERE AROUND C-ALPHA CONTAINING
1896              ALL ATOMS OF A RESIDUE
1897  CADIST   - MINIMUM DISTANCE BETWEEN ALPHA-CARBON ATOMS SUCH THAT NO
1898              BACKBONE HYDROGEN BONDS CAN BE FORMED
1899  DIST     - SMALLEST ALLOWED DISTANCE BETWEEN ANY ATOMS
1900  MAXDIST  - LARGEST ALLOWED DISTANCE BETWEEN SIDECHAIN ATOM AND C-ALPHA
1901              WITHIN A RESIDUE
1902  Q        - COUPLING CONSTANT FOR ELECTROSTATIC ENERGY
1903                     Q=-332*0.42*0.2*1000.0
1904  HBLOW    - LOWEST ALLOWED  ENERGY OF A HYDROGEN BOND IN CAL/MOL
1905  HBHIGH   - HIGHEST ALLOWED ENERGY OF A HYDROGEN BOND IN CAL/MOL
1906   --------------------------------------------------------------------*/
1907
1908 #define RN              1.65
1909 #define RCA             1.87
1910 #define RC              1.76
1911 #define RO              1.4
1912 #define RSIDEATOM       1.8
1913 #define RWATER          1.4
1914 #define SSDIST          3.0
1915 #define BREAKDIST       2.5
1916 #define RESRAD          10.0
1917 #define CADIST          9.0
1918 #define DIST            0.5
1919 #define MAXDIST         10.0
1920 #define Q               (-27888.0)
1921
1922 #define HBLOW           (-9900)
1923 #define HBHIGH          (-500)
1924
1925
1926 /***/
1927 /***************** GLOBAL DATA TYPE DEFINITIONS ***********************/
1928
1929 typedef double vector[3];
1930 typedef Char char4[4];
1931 typedef Char char6[6];
1932 typedef enum {
1933   parallel, antiparallel, nobridge
1934 } bridgetyp;
1935 typedef enum {
1936   symbol, turn3, turn4, turn5, bend, chirality, beta1, beta2
1937 } structure;
1938
1939 typedef struct hydrogenbond {
1940   long residue, energy;
1941 } hydrogenbond;
1942
1943 typedef hydrogenbond bonds[2];
1944
1945 typedef struct backbone {
1946   char6 aaident;
1947   Char sheetlabel, aa;
1948   char4 threelettercode;
1949   Char ss[(long)beta2 - (long)symbol + 1];
1950   long partner[(long)beta2 - (long)beta1 + 1];
1951   long access;
1952   double alpha, kappa;
1953   bonds acceptor, donor;
1954   vector h, n, ca, c, o;
1955   long atompointer, nsideatoms;
1956 } backbone;
1957
1958 typedef struct bridge {
1959   Char sheetname, laddername;
1960   bridgetyp btyp;
1961   long linkset[MAXBRIDGE / 32 + 2];
1962   long ib, ie, jb, je, from, towards;
1963 } bridge;
1964
1965
1966 static int bVerbose;
1967 Static int silentFlag;
1968 Static long nss, nssintra, nssinter, lchain, nbridge;
1969 Static char6 ssbonds[MAXSS][2];
1970 Static backbone chain[NMAX + 1];
1971 FILE *tapein, *tapeout;
1972 Static vector sidechain[MAXATOM];
1973 Static bridge bridgetable[MAXBRIDGE];
1974
1975 Static Void VecCopy(dest,source)
1976 double* dest;
1977 double* source;
1978 {
1979   dest[0]=source[0];
1980   dest[1]=source[1];
1981   dest[2]=source[2];
1982 }
1983
1984 Static Void StrCopy(dest,source,n)
1985 char* dest;
1986 char* source;
1987 int n;
1988 {
1989   int i;
1990   for(i=0;i<n;i++)
1991     dest[i]=source[i];
1992 }
1993
1994 Static double Atan2(y, x)
1995 double y, x;
1996 {
1997   double z;
1998
1999   if (x != 0.0)
2000     z = atan(y / x);
2001   else if (y > 0.0)
2002     z = PIHALF;
2003   else if (y < 0.0)
2004     z = -PIHALF;
2005   else
2006     z = TWOPI;
2007   if (x >= 0.0)
2008     return z;
2009   if (y > 0.0)
2010     z += PI;
2011   else
2012     z -= PI;
2013   return z;
2014 }  /* Atan2 */
2015
2016
2017 /***/
2018
2019 Static Void Diff(x, y, z)
2020 double *x, *y, *z;
2021 {
2022   z[0] = x[0] - y[0];
2023   z[1] = x[1] - y[1];
2024   z[2] = x[2] - y[2];
2025 }  /* Diff */
2026
2027
2028 /***/
2029
2030 Static double Dot(x, y)
2031 double *x, *y;
2032 {
2033   return (x[0] * y[0] + x[1] * y[1] + x[2] * y[2]);
2034 }  /* Dot */
2035
2036
2037 /***/
2038
2039 Static Void Cross(x, y, z)
2040 double *x, *y, *z;
2041 {
2042   z[0] = x[1] * y[2] - y[1] * x[2];
2043   z[1] = x[2] * y[0] - y[2] * x[0];
2044   z[2] = x[0] * y[1] - y[0] * x[1];
2045 }  /* Cross */
2046
2047
2048 /***/
2049
2050 Static Void Norm(x, xnorm)
2051 double *x;
2052 double *xnorm;
2053 {
2054   /* RETURNS INPUT VECTOR X NORMALIZED TO UNIT LENGTH.
2055      XNORM IS THE ORIGINAL LENGTH OF X.                         */
2056   double TEMP, TEMP1, TEMP2;
2057
2058   TEMP = x[0];
2059   TEMP1 = x[1];
2060   TEMP2 = x[2];
2061   *xnorm = TEMP * TEMP + TEMP1 * TEMP1 + TEMP2 * TEMP2;
2062   if (*xnorm <= 0.0)
2063     return;
2064   *xnorm = sqrt(*xnorm);
2065   x[0] /= *xnorm;
2066   x[1] /= *xnorm;
2067   x[2] /= *xnorm;
2068 }  /* Norm */
2069
2070
2071 /***/
2072
2073 Static double Dihedralangle(v1, v2, v3, v4)
2074 double *v1, *v2, *v3, *v4;
2075 {
2076   /*CALCULATES TORSION ANGLE OF A SET OF 4 ATOMS V1-V2-V3-V4.
2077     DIHEDRALANGLE IS THE ANGLE BETWEEN THE PROJECTION OF
2078     V1-V2 AND THE PROJECTION OF V4-V3 ONTO A PLANE NORMAL TO
2079     BOND V2-V3.*/
2080   /***/
2081   double Result, u, v;
2082   vector v12, v43, x, y, z, p;
2083
2084   Diff(v1, v2, v12);
2085   Diff(v4, v3, v43);
2086   Diff(v2, v3, z);
2087   Cross(z, v12, p);
2088   Cross(z, v43, x);
2089   Cross(z, x, y);
2090   u = Dot(x, x);
2091   v = Dot(y, y);
2092   Result = 360.0;
2093   if (u <= 0.0 || v <= 0.0)
2094     return Result;
2095   u = Dot(p, x) / sqrt(u);
2096   v = Dot(p, y) / sqrt(v);
2097   if (u != 0.0 || v != 0.0)
2098     return (Atan2(v, u) * RADIAN);
2099   return Result;
2100 }  /* Dihedralangle */
2101
2102
2103 /***/
2104
2105 Static double Cosangle(v1, v2, v3, v4)
2106 double *v1, *v2, *v3, *v4;
2107 {
2108   vector u, v;
2109   double x;
2110
2111   Diff(v1, v2, u);
2112   Diff(v3, v4, v);
2113   x = Dot(u, u) * Dot(v, v);
2114   if (x > 0.0)
2115     return (Dot(u, v) / sqrt(x));
2116   else
2117     return 0.0;
2118 }  /* Cosangle */
2119
2120
2121 /***/
2122
2123 Static double Distance(u, v)
2124 double *u, *v;
2125 {
2126   double TEMP, TEMP1, TEMP2;
2127
2128   TEMP = u[0] - v[0];
2129   TEMP1 = u[1] - v[1];
2130   TEMP2 = u[2] - v[2];
2131   return sqrt(TEMP * TEMP + TEMP1 * TEMP1 + TEMP2 * TEMP2);
2132 }  /* Distance */
2133
2134
2135 /***/
2136
2137 Static double Distsq(u, v)
2138 double *u, *v;
2139 {
2140   double TEMP, TEMP1, TEMP2;
2141
2142   TEMP = u[0] - v[0];
2143   TEMP1 = u[1] - v[1];
2144   TEMP2 = u[2] - v[2];
2145   return (TEMP * TEMP + TEMP1 * TEMP1 + TEMP2 * TEMP2);
2146 }  /* Distsq */
2147
2148
2149 /*--------------------------------------------------------------------*/
2150
2151 Static boolean Nochainbreak(i, j)
2152 long i, j;
2153 {
2154   long k;
2155   boolean test;
2156
2157   test = (i >= 1 && j <= NMAX && i <= j);
2158   k = i;
2159   while (test && k <= j) {
2160     if (chain[k].aa == '!')
2161       test = false;
2162     else
2163       k++;
2164   }
2165   return test;
2166 }  /* Nochainbreak */
2167
2168
2169 /***/
2170 /*--------------------------------------------------------------------*/
2171
2172 Static Void Writeresidue(res)
2173 backbone res;
2174 {
2175   long i;
2176
2177   for (i = 0; i <= 3; i++)
2178     putchar(res.threelettercode[i]);
2179   for (i = 0; i <= 5; i++)
2180     putchar(res.aaident[i]);
2181 }  /* Writeresidue */
2182
2183
2184 #define MAXSIDEATOMS    20
2185
2186
2187 typedef enum {
2188   headercard, compndcard, sourcecard, authorcard, ssbondcard, atomcard,
2189   tercard, endcard, othercard
2190 } cardtype;
2191 /***/
2192
2193 typedef struct cardcontents {
2194   cardtype art;
2195   union {
2196     Char z[128];
2197     char6 r[2];
2198     struct {
2199       char4 atomname, aaname;
2200       Char altloc, residuename;
2201       char6 reseqnum;
2202       vector coordinates;
2203     } U5;
2204     Char ch;
2205   } UU;
2206 } cardcontents;   /* CARDCONTENTS TYPE DEFINITION */
2207
2208 /***/
2209
2210
2211 Static jmp_buf _JL99;
2212
2213 /* Local variables for Inputcoordinates: */
2214 struct LOC_Inputcoordinates {
2215   long *lchain, latom, hatoms;
2216   boolean nmissing, camissing, cmissing, omissing, corelimit;
2217   vector sidecoordinates[MAXSIDEATOMS];
2218   double dco;
2219   char4 sideatomnames[MAXSIDEATOMS];
2220   backbone reszero, resinfo;
2221 } ;
2222
2223 /***/
2224
2225 Local Char Onelettercode(aaa, LINK)
2226 Char *aaa;
2227 struct LOC_Inputcoordinates *LINK;
2228 {
2229   Char aasymbol[50];
2230   Char aminoacid[150];
2231   Char string[5][30];
2232   long i, l, k;
2233   Char a;
2234
2235   StrCopy(aasymbol, "ARNDCEQGHILKMFPSTWYVBZXXXXXXXXXXXXXXXX--CCCCIPPPW-", 50L);
2236   StrCopy(string[0], "ALAARGASNASPCYSGLUGLNGLYHISILE", 30L);
2237   StrCopy(string[1], "LEULYSMETPHEPROSERTHRTRPTYRVAL", 30L);
2238   StrCopy(string[2], "ASXGLXACDALBALIABUAROBASBETHSE", 30L);
2239   StrCopy(string[3], "HYPHYLORNPCASARTAUTHYUNKACEFOR", 30L);
2240   StrCopy(string[4], "CYHCSHCSSCYXILUPRZPR0CPRTRYHOH", 30L);
2241   l = 0;
2242   for (k = 0; k <= 4; k++) {
2243     for (i = 0; i <= 29; i++) {
2244       l++;
2245       aminoacid[l - 1] = string[k][i];
2246     }
2247   }
2248   a = '-';
2249   i = 1;
2250   k = 1;
2251   while (k < 51 && a == '-') {
2252     if (aminoacid[i - 1] == aaa[0]) {
2253       if (aminoacid[i] == aaa[1]) {
2254         if (aminoacid[i + 1] == aaa[2])
2255           a = aasymbol[k - 1];
2256       }
2257     }
2258     i += 3;
2259     k++;
2260   }
2261   return a;
2262 }  /* Onelettercode */
2263
2264 /* Local variables for Checksideatoms: */
2265 struct LOC_Checksideatoms {
2266   struct LOC_Inputcoordinates *LINK;
2267 } ;
2268
2269 /***/
2270
2271 Local Void Checkdist(resinfo, LINK)
2272 backbone *resinfo;
2273 struct LOC_Checksideatoms *LINK;
2274 {
2275   long i, j, FORLIM;
2276
2277   i = 1;
2278   while (i <= resinfo->nsideatoms) {
2279     if (Distance(resinfo->ca, LINK->LINK->sidecoordinates[i - 1]) <= MAXDIST) {
2280       i++;
2281       continue;
2282     }
2283     if (bVerbose) {
2284       printf(" !!! RESIDUE ");
2285       Writeresidue(*resinfo);
2286       printf(" HAS ILLEGAL SIDECHAIN ATOM NAMED ");
2287       for (j = 0; j <= 3; j++)
2288         putchar(LINK->LINK->sideatomnames[i - 1][j]);
2289       printf(".\n");
2290       printf("     THIS ATOM WILL BE IGNORED !!!\n\n");
2291     }
2292     FORLIM = resinfo->nsideatoms;
2293     for (j = i + 1; j <= FORLIM; j++) {
2294       StrCopy(LINK->LINK->sideatomnames[j - 2],
2295               LINK->LINK->sideatomnames[j - 1], sizeof(char4));
2296       VecCopy(LINK->LINK->sidecoordinates[j - 2],
2297               LINK->LINK->sidecoordinates[j - 1]);
2298     }
2299     resinfo->nsideatoms--;
2300   }
2301 }  /* Checkdist */
2302
2303 /***/
2304
2305 Local Void Checksideatoms(resinfo, LINK)
2306 backbone *resinfo;
2307 struct LOC_Inputcoordinates *LINK;
2308 {
2309   struct LOC_Checksideatoms V;
2310   long i, j;
2311   Char c;
2312
2313   /***/
2314
2315   V.LINK = LINK;
2316   Checkdist(resinfo, &V);
2317   i = -1;
2318   c = resinfo->aa;
2319   if (c == 'G')
2320     i = 0;
2321   if (c == 'A')
2322     i = 1;
2323   if (c == 'S' || c == 'C')
2324     i = 2;
2325   if (c == 'P' || c == 'T' || c == 'V')
2326     i = 3;
2327   if (c == 'B' || c == 'M' || c == 'L' || c == 'I' || c == 'D' || c == 'N')
2328     i = 4;
2329   if (c == 'Z' || c == 'K' || c == 'Q' || c == 'E')
2330     i = 5;
2331   if (c == 'H')
2332     i = 6;
2333   if (c == 'F' || c == 'R')
2334     i = 7;
2335   if (c == 'Y')
2336     i = 8;
2337   if (c == 'W')
2338     i = 10;
2339   if ((resinfo->nsideatoms < i) && (bVerbose)) {
2340       printf(" !!! RESIDUE ");
2341       Writeresidue(*resinfo);
2342       printf(" HAS%3ld INSTEAD OF EXPECTED ", resinfo->nsideatoms);
2343       printf("%3ld SIDECHAIN ATOMS.\n", i);
2344       printf("     CALCULATED SOLVENT ACCESSIBILITY REFERS TO INCOMPLETE "
2345              "SIDECHAIN !!!\n\n");
2346     }
2347   if (i == -1 || resinfo->nsideatoms <= i)
2348     return;
2349   if (bVerbose) {
2350     printf(" !!! RESIDUE ");
2351     Writeresidue(*resinfo);
2352     printf(" HAS%3ld INSTEAD OF EXPECTED ", resinfo->nsideatoms);
2353     printf("%3ld SIDECHAIN ATOMS.\n", i);
2354     printf("     LAST SIDECHAIN ATOM NAME IS ");
2355     for (j = 0; j <= 3; j++)
2356       putchar(LINK->sideatomnames[resinfo->nsideatoms - 1][j]);
2357     printf("\n     CALCULATED SOLVENT ACCESSIBILITY INCLUDES EXTRA ATOMS !!!\n\n");
2358   }
2359 }  /* Checksideatoms */
2360
2361 /***/
2362
2363 Local Void Putresidue(LINK)
2364 struct LOC_Inputcoordinates *LINK;
2365 {
2366   /* insert residue into protein chain */
2367   long i;
2368   boolean complete;
2369   long FORLIM;
2370
2371   complete = !(LINK->nmissing || LINK->camissing || LINK->cmissing ||
2372                LINK->omissing);
2373   if (!complete &&
2374       strncmp(LINK->reszero.aaident, LINK->resinfo.aaident, sizeof(char6))
2375       && bVerbose) {
2376     printf(" !!! BACKBONE INCOMPLETE FOR RESIDUE ");
2377     Writeresidue(LINK->resinfo);
2378     printf("\n     RESIDUE WILL BE IGNORED !!!\n\n");
2379   }
2380   LINK->corelimit = (LINK->latom + LINK->resinfo.nsideatoms > MAXATOM ||
2381                      *LINK->lchain > NMAX - 2);
2382   if (complete && !LINK->corelimit) {
2383     Checksideatoms(&LINK->resinfo, LINK);
2384     VecCopy(LINK->resinfo.h, LINK->resinfo.n);
2385     if (Nochainbreak(*LINK->lchain, *LINK->lchain)) {
2386       if (Distance(chain[*LINK->lchain].c, LINK->resinfo.n) > BREAKDIST)
2387         /* keep ! at LCHAIN */
2388         /* CS Oct 1987 */
2389         if (bVerbose) {
2390           printf(" !!! EXCESSIVE C TO N DISTANCE ");
2391           printf("% .5E>% .5E\n",
2392                  Distance(chain[*LINK->lchain].c, LINK->resinfo.n), BREAKDIST);
2393           printf("     BEFORE RESIDUE ");
2394           Writeresidue(LINK->resinfo);
2395           printf(". CHAIN BREAK RESIDUE INSERTED !!!\n\n");
2396           (*LINK->lchain)++;
2397         }
2398     }
2399     if (Nochainbreak(*LINK->lchain, *LINK->lchain) && LINK->resinfo.aa != 'P') {
2400       LINK->dco = Distance(chain[*LINK->lchain].c, chain[*LINK->lchain].o);
2401       for (i = 0; i <= 2; i++)
2402         LINK->resinfo.h[i] = LINK->resinfo.n[i] +
2403             (chain[*LINK->lchain].c[i] - chain[*LINK->lchain].o[i]) / LINK->dco;
2404     }
2405     (*LINK->lchain)++;
2406     chain[*LINK->lchain] = LINK->resinfo;
2407     FORLIM = LINK->resinfo.nsideatoms;
2408     for (i = 0; i < FORLIM; i++)
2409       VecCopy(sidechain[LINK->latom + i], LINK->sidecoordinates[i]);
2410     LINK->latom += LINK->resinfo.nsideatoms;
2411   }
2412   if (Nochainbreak(*LINK->lchain, *LINK->lchain) && !complete)
2413     (*LINK->lchain)++;
2414   LINK->resinfo = LINK->reszero;
2415   LINK->nmissing = true;
2416   LINK->camissing = true;
2417   LINK->cmissing = true;
2418   LINK->omissing = true;
2419 }  /* Putresidue */
2420
2421 /***/
2422
2423 Local Void Getresidue(atomname, coordinates, LINK)
2424 Char *atomname;
2425 double *coordinates;
2426 struct LOC_Inputcoordinates *LINK;
2427 {
2428   boolean hydrogenatom;
2429
2430   hydrogenatom = ((atomname[0] == '9' || atomname[0] == '8' ||
2431                    atomname[0] == '7' || atomname[0] == '6' ||
2432                    atomname[0] == '5' || atomname[0] == '4' ||
2433                    atomname[0] == '3' || atomname[0] == '2' ||
2434                    atomname[0] == '1' || atomname[0] == '0' ||
2435                    atomname[0] == ' ') &&
2436                   (atomname[1] == 'D' || atomname[1] == 'H'));
2437   if (hydrogenatom) {
2438     LINK->hatoms++;
2439     return;
2440   }
2441   if (!strncmp(atomname, " N  ", sizeof(char4))) {
2442     LINK->nmissing = false;
2443     VecCopy(LINK->resinfo.n, coordinates);
2444     return;
2445   }
2446   if (!strncmp(atomname, " CA ", sizeof(char4))) {
2447     LINK->camissing = false;
2448     VecCopy(LINK->resinfo.ca, coordinates);
2449     return;
2450   }
2451   if (!strncmp(atomname, " C  ", sizeof(char4))) {
2452     LINK->cmissing = false;
2453     VecCopy(LINK->resinfo.c, coordinates);
2454     return;
2455   }
2456   if (!strncmp(atomname, " O  ", sizeof(char4))) {
2457     LINK->omissing = false;
2458     VecCopy(LINK->resinfo.o, coordinates);
2459     return;
2460   }
2461   if (LINK->resinfo.nsideatoms >= MAXSIDEATOMS)
2462     return;
2463   LINK->resinfo.nsideatoms++;
2464   VecCopy(LINK->sidecoordinates[LINK->resinfo.nsideatoms - 1], coordinates);
2465   StrCopy(LINK->sideatomnames[LINK->resinfo.nsideatoms - 1], atomname,
2466          sizeof(char4));
2467 }  /* Getresidue */
2468
2469 /***/
2470
2471 Local Void Readcard(cardinfo, LINK)
2472 cardcontents *cardinfo;
2473 struct LOC_Inputcoordinates *LINK;
2474 {
2475   Char c;
2476   long k, l, m;
2477   char6 key;
2478
2479   cardinfo->art = othercard;
2480   do {
2481     if (!P_eof(tapein)) {
2482       *key = getc(tapein);
2483       if (key[0] == '\n')
2484         key[0] = ' ';
2485     }
2486   } while (!(isupper(key[0]) | P_eof(tapein)));
2487   if (P_eof(tapein)) {
2488     cardinfo->art = endcard;
2489     return;
2490   }
2491   for (l = 1; l <= 5; l++) {
2492     if (!P_eoln(tapein)) {
2493       key[l] = getc(tapein);
2494       if (key[l] == '\n')
2495         key[l] = ' ';
2496     }
2497   }
2498   if (!strncmp(key, "HEADER", sizeof(char6)))
2499     cardinfo->art = headercard;
2500   if (!strncmp(key, "COMPND", sizeof(char6)))
2501     cardinfo->art = compndcard;
2502   if (!strncmp(key, "SOURCE", sizeof(char6)))
2503     cardinfo->art = sourcecard;
2504   if (!strncmp(key, "AUTHOR", sizeof(char6)))
2505     cardinfo->art = authorcard;
2506   if (!strncmp(key, "SSBOND", sizeof(char6)))
2507     cardinfo->art = ssbondcard;
2508   if (!strncmp(key, "ATOM  ", sizeof(char6)))
2509     cardinfo->art = atomcard;
2510   if (!strncmp(key, "TER   ", sizeof(char6)))
2511     cardinfo->art = tercard;
2512   if (!strncmp(key, "END   ", sizeof(char6)))
2513     cardinfo->art = endcard;
2514   switch (cardinfo->art) {
2515
2516   case headercard:
2517   case compndcard:
2518   case sourcecard:
2519   case authorcard:
2520     for (l = 0; l <= 5; l++)
2521       cardinfo->UU.z[l] = key[l];
2522     for (l = 6; l <= 126; l++)
2523       cardinfo->UU.z[l] = ' ';
2524     cardinfo->UU.z[127] = '.';
2525     if (cardinfo->art == headercard)
2526       m = 66;
2527     else
2528       m = 70;
2529     for (l = 6; l < m; l++) {
2530       if (!P_eoln(tapein)) {
2531         cardinfo->UU.z[l] = getc(tapein);
2532         if (cardinfo->UU.z[l] == '\n')
2533           cardinfo->UU.z[l] = ' ';
2534       }
2535     }
2536     break;
2537
2538   case ssbondcard:
2539     for (l = 7; l <= 8; l++) {
2540       c = getc(tapein);
2541       if (c == '\n')
2542         c = ' ';
2543     }
2544     for (k = 0; k <= 1; k++) {
2545       for (l = 1; l <= 7; l++) {
2546         c = getc(tapein);
2547         if (c == '\n')
2548           c = ' ';
2549       }
2550       cardinfo->UU.r[k][5] = getc(tapein);
2551       c = getc(tapein);
2552       if (cardinfo->UU.r[k][5] == '\n')
2553         cardinfo->UU.r[k][5] = ' ';
2554       if (c == '\n')
2555         c = ' ';
2556       /* minor modification suggested by Steven Sheriff */
2557       for (l = 0; l <= 3; l++) {
2558         cardinfo->UU.r[k][l] = getc(tapein);
2559         if (cardinfo->UU.r[k][l] == '\n')
2560           cardinfo->UU.r[k][l] = ' ';
2561       }
2562       if (P_eoln(tapein))
2563         cardinfo->UU.r[k][4] = ' ';
2564       else {
2565         cardinfo->UU.r[k][4] = getc(tapein);
2566         if (cardinfo->UU.r[k][4] == '\n')
2567           cardinfo->UU.r[k][4] = ' ';
2568       }
2569     }
2570     /* end minor modification suggested by Steven Sheriff */
2571     break;
2572
2573   case atomcard:
2574     for (l = 7; l <= 12; l++) {
2575       c = getc(tapein);
2576       if (c == '\n')
2577         c = ' ';
2578     }
2579     for (l = 0; l <= 3; l++) {
2580       cardinfo->UU.U5.atomname[l] = getc(tapein);
2581       if (cardinfo->UU.U5.atomname[l] == '\n')
2582         cardinfo->UU.U5.atomname[l] = ' ';
2583     }
2584     cardinfo->UU.U5.altloc = getc(tapein);
2585     if (cardinfo->UU.U5.altloc == '\n')
2586       cardinfo->UU.U5.altloc = ' ';
2587     for (l = 0; l <= 2; l++) {
2588       cardinfo->UU.U5.aaname[l] = getc(tapein);
2589       if (cardinfo->UU.U5.aaname[l] == '\n')
2590         cardinfo->UU.U5.aaname[l] = ' ';
2591     }
2592     cardinfo->UU.U5.aaname[3] = ' ';
2593     cardinfo->UU.U5.residuename = Onelettercode(cardinfo->UU.U5.aaname, LINK);
2594     c = getc(tapein);
2595     cardinfo->UU.U5.reseqnum[5] = getc(tapein);
2596     if (c == '\n')
2597       c = ' ';
2598     if (cardinfo->UU.U5.reseqnum[5] == '\n')
2599       cardinfo->UU.U5.reseqnum[5] = ' ';
2600     for (l = 0; l <= 4; l++) {
2601       cardinfo->UU.U5.reseqnum[l] = getc(tapein);
2602       if (cardinfo->UU.U5.reseqnum[l] == '\n')
2603         cardinfo->UU.U5.reseqnum[l] = ' ';
2604     }
2605     for (l = 0; l <= 2; l++)
2606       fscanf(tapein, "%lf", &cardinfo->UU.U5.coordinates[l]);
2607     break;
2608
2609   case tercard:
2610   case endcard:
2611   case othercard:
2612     /* blank case */
2613     break;
2614   }
2615   fscanf(tapein, "%*[^\n]");
2616   getc(tapein);
2617 }  /* Readcard */
2618
2619
2620 /***/
2621 /*--------------------------------------------------------------------*/
2622 /* SEE BROOKHAVEN PROTEIN DATA BANK ATOMIC COORDINATE ENTRY FORMAT
2623                                     OF DEC. 1981.
2624    -------------------------------------------------------------------*/
2625
2626 Static Void Inputcoordinates(lchain_)
2627 long *lchain_;
2628 {
2629   struct LOC_Inputcoordinates V;
2630   Char datestring[30];
2631   long i, j;
2632   boolean finish;
2633   structure s;
2634   cardtype ctype;
2635   cardcontents cardinfo;
2636   long cardhist[(long)othercard - (long)headercard + 1];
2637
2638   /***/
2639
2640   V.lchain = lchain_;
2641   nss = 0;
2642   V.latom = 0;
2643   V.hatoms = 0;
2644   for (j = 0; j <= 5; j++)
2645     V.reszero.aaident[j] = ' ';
2646   V.reszero.aa = '!';
2647   V.reszero.access = 0;
2648   StrCopy(V.reszero.threelettercode, "    ", sizeof(char4));
2649   for (s = symbol; (long)s <= (long)beta2; s = (structure)((long)s + 1))
2650     V.reszero.ss[(long)s - (long)symbol] = ' ';
2651   V.reszero.sheetlabel = ' ';
2652   V.reszero.partner[0] = 0;
2653   V.reszero.partner[(long)beta2 - (long)beta1] = 0;
2654   V.reszero.alpha = 360.0;
2655   V.reszero.kappa = 360.0;
2656   for (j = 0; j <= 1; j++) {
2657     V.reszero.acceptor[j].residue = 0;
2658     V.reszero.acceptor[j].energy = 0;
2659     V.reszero.donor[j].residue = 0;
2660     V.reszero.donor[j].energy = 0;
2661   }
2662   V.reszero.atompointer = 0;
2663   V.reszero.nsideatoms = 0;
2664   for (j = 0; j <= 2; j++) {
2665     V.reszero.h[j] = 0.0;
2666     V.reszero.n[j] = 0.0;
2667     V.reszero.ca[j] = 0.0;
2668     V.reszero.c[j] = 0.0;
2669     V.reszero.o[j] = 0.0;
2670   }
2671   for (i = 0; i <= NMAX; i++)
2672     chain[i] = V.reszero;
2673   Date(datestring);   /* DATE(DAY-MONTH-YEAR); */
2674   /* comment out this line if necessary */
2675   fprintf(tapeout, "**** SECONDARY STRUCTURE DEFINITION ");
2676   fprintf(tapeout, "BY THE PROGRAM DSSP, VERSION OCT. 1988 ****");
2677   fprintf(tapeout, " DATE=%.11s", datestring);
2678   for (i = 106; i <= 127; i++)
2679     putc(' ', tapeout);
2680   fprintf(tapeout, ".\n");
2681   fprintf(tapeout, "REFERENCE W. KABSCH AND C.SANDER, BIOPOLYMERS ");
2682   fprintf(tapeout, "22 (1983) 2577-2637");
2683   for (i = 66; i <= 127; i++)
2684     putc(' ', tapeout);
2685   fprintf(tapeout, ".\n");
2686   for (ctype = headercard;
2687        (long)ctype <= (long)othercard;
2688        ctype = (cardtype)((long)ctype + 1))
2689     cardhist[(long)ctype - (long)headercard] = 0;
2690   V.corelimit = false;
2691   finish = false;
2692   V.resinfo = V.reszero;
2693   V.nmissing = true;
2694   V.camissing = true;
2695   V.cmissing = true;
2696   V.omissing = true;
2697   do {
2698     Readcard(&cardinfo, &V);
2699     cardhist[(long)cardinfo.art - (long)headercard]++;
2700     switch (cardinfo.art) {
2701
2702     case headercard:
2703     case compndcard:
2704     case sourcecard:
2705     case authorcard:
2706       if (cardhist[(long)cardinfo.art - (long)headercard] == 1) {
2707         for (i = 0; i <= 127; i++)
2708           putc(cardinfo.UU.z[i], tapeout);
2709         putc('\n', tapeout);
2710       }
2711       break;
2712
2713     case ssbondcard:
2714       nss++;
2715       for (i = 0; i <= 1; i++)
2716         StrCopy(ssbonds[nss - 1][i], cardinfo.UU.r[i], sizeof(char6));
2717       break;
2718
2719     case atomcard:
2720       if (cardinfo.UU.U5.residuename != '-' &&
2721           (cardinfo.UU.U5.altloc == 'A' || cardinfo.UU.U5.altloc == ' ')) {
2722         if (strncmp(V.resinfo.aaident, cardinfo.UU.U5.reseqnum, sizeof(char6))) {
2723           Putresidue(&V);
2724           V.resinfo.atompointer = V.latom;
2725           StrCopy(V.resinfo.aaident, cardinfo.UU.U5.reseqnum, sizeof(char6));
2726           V.resinfo.aa = cardinfo.UU.U5.residuename;
2727           StrCopy(V.resinfo.threelettercode, cardinfo.UU.U5.aaname,
2728                  sizeof(char4));
2729         }
2730         Getresidue(cardinfo.UU.U5.atomname, cardinfo.UU.U5.coordinates, &V);
2731       }
2732       if ((cardinfo.UU.U5.residuename == '-') && bVerbose) {
2733         printf(" !!! RESIDUE ");
2734         for (i = 0; i <= 3; i++)
2735           putchar(cardinfo.UU.U5.aaname[i]);
2736         for (i = 0; i <= 5; i++)
2737           putchar(cardinfo.UU.U5.reseqnum[i]);
2738         printf(" HAS NONSTANDARD NAME.\n");
2739         printf("     RESIDUE WILL BE ");
2740         printf("IGNORED !!!\n");
2741       }
2742       if ((cardinfo.UU.U5.altloc != 'A' && cardinfo.UU.U5.altloc != ' ') &&
2743           bVerbose) {
2744         printf(" !!! IN RESIDUE");
2745         for (i = 0; i <= 3; i++)
2746           printf(" %c", cardinfo.UU.U5.aaname[i]);
2747         for (i = 0; i <= 5; i++)
2748           putchar(cardinfo.UU.U5.reseqnum[i]);
2749         printf(" ALTERNATE LOCATION INDICATOR ");
2750         printf("IS %c AND\n", cardinfo.UU.U5.altloc);
2751         printf("     NOT BLANK OR A. ATOM ");
2752         printf("NAMED ");
2753         for (i = 0; i <= 3; i++)
2754           putchar(cardinfo.UU.U5.atomname[i]);
2755         printf(" WILL BE IGNORED !!!\n\n");
2756       }
2757       break;
2758
2759     case tercard:
2760       Putresidue(&V);
2761       break;
2762
2763     case endcard:
2764       finish = true;
2765       Putresidue(&V);
2766       break;
2767
2768     case othercard:
2769       /* blank case */
2770       break;
2771     }
2772   } while (!(V.corelimit || finish));
2773   if ((V.corelimit) && bVerbose) {
2774     printf(" !!! NUMBER OF ATOMS OR RESIDUES EXCEEDS ");
2775     printf("STORAGE CAPACITY !!!\n");
2776   }
2777   if (!Nochainbreak(*V.lchain, *V.lchain))
2778     (*V.lchain)--;
2779   if ((V.hatoms > 0) && bVerbose) {
2780     printf(" !!! %12ld HYDROGEN OR DEUTERIUM ATOMS WERE IGNORED\n", V.hatoms);
2781     printf("     IN THE CALCULATION OF SIDE CHAIN SOLVENT \n");
2782     printf("     ACCESSIBILITY !!!\n");
2783   }
2784   if (bVerbose) {
2785     if (cardhist[0] < 1)
2786       printf(" !!! HEADER-CARD MISSING !!!\n");
2787     if (cardhist[(long)compndcard - (long)headercard] < 1)
2788       printf(" !!! COMPOUND-CARD MISSING !!!\n");
2789     if (cardhist[(long)sourcecard - (long)headercard] < 1)
2790       printf(" !!! SOURCE-CARD MISSING !!!\n");
2791     if (cardhist[(long)authorcard - (long)headercard] < 1)
2792       printf(" !!! AUTHOR CARD MISSING !!!\n");
2793     if (*V.lchain < 1) {
2794       printf(" !!! NO RESIDUE WITH COMPLETE BACKBONE !!!\n");
2795       longjmp(_JL99, 1);
2796     }
2797     if (V.latom == 0)
2798       printf(" !!! ALL SIDECHAIN COORDINATES MISSING !!!\n");
2799   }
2800 }  /* Inputcoordinates */
2801
2802 #undef MAXSIDEATOMS
2803
2804
2805 /***/
2806 /*--------------------------------------------------------------------*/
2807
2808 Static boolean Testbond(i, j)
2809 long i, j;
2810 {
2811   /* TESTBOND IS TRUE IF I IS DONOR[=NH] TO J, OTHERWISE FALSE */
2812   backbone *WITH;
2813
2814   WITH = &chain[i];
2815   return (WITH->acceptor[0].residue == j && WITH->acceptor[0].energy < HBHIGH ||
2816           WITH->acceptor[1].residue == j && WITH->acceptor[1].energy < HBHIGH);
2817 }  /* Testbond */
2818
2819
2820 /***/
2821
2822 Local boolean Testssbond(i, j)
2823 long i, j;
2824 {
2825   boolean ssbond;
2826   long k;
2827
2828   ssbond = false;
2829   k = 1;
2830   if (!(Nochainbreak(i, i) & Nochainbreak(j, j)))
2831     return ssbond;
2832   while (!(ssbond || k > nss)) {
2833     ssbond = ((!strncmp(chain[i].aaident, ssbonds[k - 1][0], sizeof(char6)) &&
2834                !strncmp(chain[j].aaident, ssbonds[k - 1][1], sizeof(char6))) ||
2835               (!strncmp(chain[i].aaident, ssbonds[k - 1][1], sizeof(char6)) &&
2836                !strncmp(chain[j].aaident, ssbonds[k - 1][0], sizeof(char6))));
2837     k++;
2838   }
2839   return ssbond;
2840 }  /* Testssbond */
2841
2842
2843 /***/
2844 /*--------------------------------------------------------------------*/
2845
2846 Static Void Flagssbonds()
2847 {
2848   boolean ssbond;
2849   Char cc;
2850   long i, j, ii, jj;
2851   double d;
2852   long FORLIM;
2853   backbone *WITH;
2854   long FORLIM1;
2855
2856   /***/
2857
2858   nssintra = 0;
2859   nssinter = 0;
2860   cc = '`';
2861   FORLIM = lchain - 2;
2862   for (i = 1; i <= FORLIM; i++) {
2863     if (chain[i].aa == 'C' && chain[i].nsideatoms > 1) {
2864       ii = chain[i].atompointer + 2;
2865       j = i + 1;
2866       do {
2867         j++;
2868         ssbond = false;
2869         if (chain[j].nsideatoms > 1 && chain[j].aa == 'C')
2870           jj = chain[j].atompointer + 2;
2871         else
2872           jj = 0;
2873         if (jj > 0)
2874           ssbond = (Distance(sidechain[ii - 1], sidechain[jj - 1]) < SSDIST);
2875       } while (!(ssbond || j == lchain));
2876       if (ssbond & (!Testssbond(i, j)))
2877         if (bVerbose) {
2878           printf(" !!! ADDITIONAL SSBOND FOUND BETWEEN ");
2879           printf("RESIDUES ");
2880           Writeresidue(chain[i]);
2881           printf(" AND ");
2882           Writeresidue(chain[j]);
2883           printf(" !!!\n\n");
2884         }
2885     }
2886   }
2887   if (nss > 0) {
2888     FORLIM = lchain - 2;
2889     for (i = 1; i <= FORLIM; i++) {
2890       WITH = &chain[i];
2891       if (WITH->aa == 'C') {
2892         FORLIM1 = lchain;
2893         for (j = i + 2; j <= FORLIM1; j++) {
2894           if (chain[j].aa == 'C') {
2895             if (Testssbond(i, j)) {
2896               if (cc == 'z') {
2897                 if (bVerbose)
2898                   printf(" !!! SS-BRIDGE LABEL RESTART AT a !!!\n");
2899                 cc = '`';
2900               }
2901               cc++;
2902               WITH->aa = cc;
2903               chain[j].aa = cc;
2904               if (Nochainbreak(i, j))
2905                 nssintra++;
2906               else
2907                 nssinter++;
2908               if (WITH->nsideatoms > 1) {
2909                 if (chain[j].nsideatoms > 1) {
2910                   jj = chain[j].atompointer + 2;
2911                   ii = WITH->atompointer + 2;
2912                   d = Distance(sidechain[ii - 1], sidechain[jj - 1]);
2913                   if ((d > SSDIST) && bVerbose) {
2914                     printf(" !!! SSBOND DISTANCE IS%5.1f BETWEEN RESIDUES", d);
2915                     Writeresidue(chain[i]);
2916                     printf(" AND ");
2917                     Writeresidue(chain[j]);
2918                     printf(" !!!\n\n");
2919                   }
2920                 }
2921               }
2922             }
2923           }
2924         }
2925       }
2926     }
2927   }
2928   if ((nss != nssintra + nssinter) && bVerbose)
2929     printf(" !!! ERROR IN SSBOND DATA RECORDS !!!\n");
2930 }  /* Flagssbonds */
2931
2932
2933 /***/
2934 /*--------------------------------------------------------------------*/
2935
2936 Static Void Flagchirality()
2937 {
2938   long i;
2939   double ckap, skap;
2940   long FORLIM;
2941   backbone *WITH;
2942
2943   FORLIM = lchain - 2;
2944   for (i = 2; i <= FORLIM; i++) {
2945     WITH = &chain[i];
2946     if (Nochainbreak(i - 1, i + 2)) {
2947       WITH->alpha = Dihedralangle(chain[i - 1].ca, WITH->ca, chain[i + 1].ca,
2948                                   chain[i + 2].ca);
2949       if (WITH->alpha < 0.0)
2950         WITH->ss[(long)chirality - (long)symbol] = '-';
2951       else
2952         WITH->ss[(long)chirality - (long)symbol] = '+';
2953     }
2954   }
2955   FORLIM = lchain - 2;
2956   /***/
2957   for (i = 3; i <= FORLIM; i++) {
2958     WITH = &chain[i];
2959     if (Nochainbreak(i - 2, i + 2)) {
2960       ckap = Cosangle(chain[i].ca, chain[i - 2].ca, chain[i + 2].ca,
2961                       chain[i].ca);
2962       skap = sqrt(1.0 - ckap * ckap);
2963       WITH->kappa = RADIAN * Atan2(skap, ckap);
2964     }
2965   }
2966 }  /* Flagchirality */
2967
2968
2969 /***/
2970
2971 Local long Bondenergy(i, j)
2972 long i, j;
2973 {
2974   /*RESIDUE I IS DONOR[=NH],J IS ACCEPTOR[=CO] OF THE PROTON IN THE
2975      HYDROGEN BOND. THE BONDENERGY IS IN CAL/MOL */
2976   double dho, dhc, dnc, dno;
2977   long hbe;
2978   backbone *WITH;
2979
2980   hbe = 0;
2981   WITH = &chain[i];
2982   if (WITH->aa == 'P')
2983     return hbe;
2984   dho = Distance(WITH->h, chain[j].o);
2985   dhc = Distance(WITH->h, chain[j].c);
2986   dnc = Distance(WITH->n, chain[j].c);
2987   dno = Distance(WITH->n, chain[j].o);
2988   if (dho < DIST || dhc < DIST || dnc < DIST || dno < DIST)
2989     hbe = HBLOW;
2990   else
2991     hbe = (long)floor(Q / dho - Q / dhc + Q / dnc - Q / dno + 0.5);
2992   if (hbe > HBLOW)
2993     return hbe;
2994   if (bVerbose) {
2995     printf(" !!! CONTACT BETWEEN RESIDUES ");
2996     Writeresidue(chain[i]);
2997     printf(" AND ");
2998     Writeresidue(chain[j]);
2999     printf("  TOO CLOSE !!!\n");
3000   }
3001   hbe = HBLOW;
3002   return hbe;
3003 }  /* Bondenergy */
3004
3005 /***/
3006
3007 Local Void Updatebonds(b, hb)
3008 hydrogenbond *b;
3009 hydrogenbond hb;
3010 {
3011   if (hb.energy < b[0].energy) {
3012     b[1] = b[0];
3013     b[0] = hb;
3014   } else if (hb.energy < b[1].energy)
3015     b[1] = hb;
3016 }  /* Updatebonds */
3017
3018 /***/
3019
3020 Local Void Setbonds(i, j)
3021 long i, j;
3022 {
3023   /*I IS NH, J IS CO*/
3024   hydrogenbond hb;
3025
3026   hb.energy = Bondenergy(i, j);
3027   hb.residue = j;
3028   /* CO(J) IS ACCEPTOR OF NH(I) */
3029   Updatebonds(chain[i].acceptor, hb);
3030   hb.residue = i;
3031   Updatebonds(chain[j].donor, hb);
3032 }  /* Setbond */
3033
3034
3035 /***/
3036 /*--------------------------------------------------------------------*/
3037
3038 Static Void Flaghydrogenbonds()
3039 {
3040   long i, j, FORLIM;
3041   backbone *WITH;
3042   long FORLIM1;
3043
3044   /***/
3045
3046   FORLIM = lchain;
3047   for (i = 1; i <= FORLIM; i++) {
3048     if (Nochainbreak(i, i)) {
3049       WITH = &chain[i];
3050       FORLIM1 = lchain;
3051       for (j = i + 1; j <= FORLIM1; j++) {
3052         if (Nochainbreak(j, j)) {
3053           if (Distance(WITH->ca, chain[j].ca) < CADIST) {
3054             Setbonds(i, j);
3055             if (j != i + 1)
3056               Setbonds(j, i);
3057           }
3058         }
3059       }
3060     }
3061   }
3062 }  /* Flaghydrogenbonds */
3063
3064
3065 /***/
3066
3067 Local Void Ladder(i, j, b)
3068 long i, j;
3069 bridgetyp b;
3070 {
3071   long k;
3072   boolean found;
3073   bridge *WITH;
3074
3075   found = false;
3076   k = 1;
3077   if (b == nobridge || i >= j)
3078     return;
3079   do {
3080     WITH = &bridgetable[k - 1];
3081     if (WITH->ib == 0) {
3082       WITH->ib = i;
3083       WITH->ie = i;
3084       WITH->jb = j;
3085       WITH->je = j;
3086       WITH->from = 0;
3087       WITH->towards = 0;
3088       WITH->btyp = b;
3089       nbridge++;
3090       found = true;
3091     } else {
3092       found = (WITH->btyp == b && i == WITH->ie + 1) & Nochainbreak(WITH->ie,
3093                 i) & (((j == WITH->je + 1 && b == parallel) &
3094                        Nochainbreak(WITH->je, j)) | ((j == WITH->jb - 1 &&
3095                           b == antiparallel) & Nochainbreak(j, WITH->jb)));
3096 /* p2c: dssp.p, line 1609: Note:
3097  * Line breaker spent 1.1+0.26 seconds, 3126 tries on line 1540 [251] */
3098       if (found) {
3099         WITH->ie++;
3100         if (b == parallel)
3101           WITH->je++;
3102         else
3103           WITH->jb--;
3104       } else {
3105         k++;
3106         if (k > MAXBRIDGE) {
3107           if (bVerbose)
3108             printf(" !!! BRIDGETABLE OVERFLOW !!!\n");
3109           longjmp(_JL99, 1);
3110         }
3111       }
3112     }
3113   } while (!found);   /* Ladder */
3114 }
3115
3116 /***/
3117
3118 Local Void Testbridge(i)
3119 long i;
3120 {
3121   long j1, j2, j;
3122   bridgetyp b;
3123
3124   /***/
3125
3126   j1 = 0;
3127   j2 = 0;
3128   j = i + 3;
3129   if (!Nochainbreak(i - 1, i + 1))
3130     return;
3131   while (j2 == 0 && j < lchain) {
3132     if (Nochainbreak(j - 1, j + 1)) {
3133       if ((Testbond(i + 1, j) & Testbond(j, i - 1)) |
3134           (Testbond(j + 1, i) & Testbond(i, j - 1)))
3135         b = parallel;
3136       else if ((Testbond(i + 1, j - 1) & Testbond(j + 1, i - 1)) |
3137                (Testbond(j, i) & Testbond(i, j)))
3138         b = antiparallel;
3139       else
3140         b = nobridge;
3141       if (b != nobridge) {
3142         if (j1 == 0) {
3143           j1 = j;
3144           Ladder(i, j, b);
3145         } else if (j != j1) {
3146           j2 = j;
3147           Ladder(i, j, b);
3148         }
3149       }
3150     }
3151     j++;
3152   }
3153 }  /* Testbridge */
3154
3155 /***/
3156
3157 Local Void Extendladder()
3158 {
3159   long i, j, ib1, jb1, je1;
3160   boolean bulge;
3161   long FORLIM;
3162   bridge *WITH;
3163   long SET[11];
3164
3165   FORLIM = nbridge;
3166   for (i = 1; i <= FORLIM; i++) {
3167     WITH = &bridgetable[i - 1];
3168     j = i + 1;
3169     while (j <= nbridge && WITH->towards == 0) {
3170       ib1 = bridgetable[j - 1].ib;
3171       jb1 = bridgetable[j - 1].jb;
3172       je1 = bridgetable[j - 1].je;
3173       bulge = (Nochainbreak(WITH->ie, ib1) && ib1 - WITH->ie < 6 &&
3174                bridgetable[j - 1].btyp == WITH->btyp &&
3175                bridgetable[j - 1].from == 0);
3176       if (bulge) {
3177         switch (WITH->btyp) {
3178
3179         case parallel:
3180           bulge = (jb1 - WITH->je < 6 && ib1 - WITH->ie < 3 ||
3181                    jb1 - WITH->je < 3) & Nochainbreak(WITH->je, jb1);
3182           break;
3183
3184         case antiparallel:
3185           bulge = (WITH->jb - je1 < 6 && ib1 - WITH->ie < 3 ||
3186                    WITH->jb - je1 < 3) & Nochainbreak(je1, WITH->jb);
3187           break;
3188         }
3189       }
3190       if (bulge) {
3191         WITH->towards = j;
3192         bridgetable[j - 1].from = i;
3193       }
3194       j++;
3195     }
3196   }
3197   FORLIM = nbridge;
3198   for (i = 1; i <= FORLIM; i++) {
3199     WITH = &bridgetable[i - 1];
3200     if (WITH->from == 0) {
3201       P_expset(WITH->linkset, 0L);
3202       j = i;
3203       do {
3204         P_addset(WITH->linkset, (int)j);
3205         j = bridgetable[j - 1].towards;
3206       } while (j != 0);
3207       j = WITH->towards;
3208       while (j != 0) {
3209         P_setcpy(bridgetable[j - 1].linkset, WITH->linkset);
3210         j = bridgetable[j - 1].towards;
3211       }
3212     }
3213   }
3214 }  /* Extendladder */
3215
3216 /* Local variables for Sheet: */
3217 struct LOC_Sheet {
3218   long ladderset[MAXBRIDGE / 32 + 2], sheetset[MAXBRIDGE / 32 + 2];
3219 } ;
3220
3221 /***/
3222
3223 Local boolean Link(l1, l2)
3224 long l1, l2;
3225 {
3226   /* LINK IS TRUE IF THERE IS A COMMON RESIDUE IN LADDERS L1 AND L2 */
3227   long ib1, ie1, jb1, je1, ib2, ie2, jb2, je2;
3228
3229   ib1 = bridgetable[l1 - 1].ib;
3230   ie1 = bridgetable[l1 - 1].ie;
3231   jb1 = bridgetable[l1 - 1].jb;
3232   je1 = bridgetable[l1 - 1].je;
3233   ib2 = bridgetable[l2 - 1].ib;
3234   ie2 = bridgetable[l2 - 1].ie;
3235   jb2 = bridgetable[l2 - 1].jb;
3236   je2 = bridgetable[l2 - 1].je;
3237   return (ie1 >= ib2 && ib1 <= ie2 || ie1 >= jb2 && ib1 <= je2 ||
3238           je1 >= ib2 && jb1 <= ie2 || je1 >= jb2 && jb1 <= je2);
3239 }  /* Link */
3240
3241 /***/
3242
3243 Local Void Findsheet(LINK)
3244 struct LOC_Sheet *LINK;
3245 {
3246   long l1, l2;
3247   boolean finish;
3248   long FORLIM, FORLIM1;
3249
3250   /***/
3251
3252   P_expset(LINK->sheetset, 0L);
3253   l1 = 0;
3254   if (*LINK->ladderset != 0L) {
3255     do {
3256       l1++;
3257     } while (!P_inset((int)l1, LINK->ladderset));
3258   }
3259   if (l1 > 0)
3260     P_setcpy(LINK->sheetset, bridgetable[l1 - 1].linkset);
3261   if (l1 <= 0)
3262     return;
3263   do {
3264     finish = true;
3265     FORLIM = nbridge;
3266     for (l1 = 1; l1 <= FORLIM; l1++) {
3267       if (P_inset((int)l1, LINK->sheetset)) {
3268         FORLIM1 = nbridge;
3269         for (l2 = 1; l2 <= FORLIM1; l2++) {
3270           if (P_inset((int)l2, LINK->ladderset)) {
3271             if (Link(l1, l2)) {
3272               P_setunion(LINK->sheetset, LINK->sheetset,
3273                          bridgetable[l2 - 1].linkset);
3274               P_setdiff(LINK->ladderset, LINK->ladderset,
3275                         bridgetable[l2 - 1].linkset);
3276               finish = false;
3277             }
3278           }
3279         }
3280       }
3281     }
3282   } while (!finish);   /* Findsheet */
3283 }
3284
3285 /***/
3286
3287 Local Void Sheet()
3288 {
3289   struct LOC_Sheet V;
3290   long asci, i, j;
3291   Char ccs;
3292   long SET[11];
3293   long FORLIM;
3294   bridge *WITH;
3295
3296   /***/
3297
3298   P_expset(V.ladderset, 0L);
3299   FORLIM = nbridge;
3300   for (i = 1; i <= FORLIM; i++)
3301     P_addset(V.ladderset, (int)i);
3302   ccs = '@';
3303   asci = 64;
3304   while (*V.ladderset != 0L) {
3305     ccs++;
3306     if (ccs > 'z') {
3307       if (bVerbose)
3308         printf(" !!! SHEET LABEL RESTART AT A !!!\n");
3309       ccs = 'A';
3310     }
3311     Findsheet(&V);
3312     FORLIM = nbridge;
3313     for (i = 1; i <= FORLIM; i++) {
3314       WITH = &bridgetable[i - 1];
3315       if (P_inset((int)i, V.sheetset) && WITH->from == 0) {
3316         if (asci == 90) {
3317           if (bVerbose)
3318             printf(" !!! STRAND LABEL RESTART AT A !!!\n");
3319           asci = 64;
3320         }
3321         asci++;
3322         if (WITH->btyp == parallel)
3323           WITH->laddername = (Char)(asci + 32);
3324         else
3325           WITH->laddername = (Char)asci;
3326         WITH->sheetname = ccs;
3327         P_setcpy(WITH->linkset, V.sheetset);
3328         j = WITH->towards;
3329         while (j != 0) {
3330           bridgetable[j - 1].laddername = WITH->laddername;
3331           bridgetable[j - 1].sheetname = WITH->sheetname;
3332           P_setcpy(bridgetable[j - 1].linkset, V.sheetset);
3333           j = bridgetable[j - 1].towards;
3334         }
3335       }
3336     }
3337   }
3338 }  /* Sheet */
3339
3340 /***/
3341
3342 Local Void Markstrands()
3343 {
3344   long i, j, l, ib0, ie0, jb0, je0;
3345   structure beta, betai, betaj;
3346   long iset[(long)beta2 - (long)beta1 + 1][9],
3347        jset[(long)beta2 - (long)beta1 + 1][9];
3348   Char cc;
3349   long FORLIM, FORLIM1;
3350   long SET[9];
3351   bridge *WITH;
3352   backbone *WITH1;
3353   long SET1[9];
3354   long SET2[3];
3355   long SET3[3];
3356
3357   FORLIM = nbridge;
3358   for (i = 1; i <= FORLIM; i++) {
3359     if (bridgetable[i - 1].from == 0) {
3360       j = i;
3361       for (beta = beta1;
3362            (long)beta <= (long)beta2;
3363            beta = (structure)((long)beta + 1)) {
3364         P_setcpy(iset[(long)beta - (long)beta1], P_expset(SET, 0L));
3365         P_setcpy(jset[(long)beta - (long)beta1], P_expset(SET, 0L));
3366       }
3367       ib0 = lchain;
3368       ie0 = 0;
3369       jb0 = lchain;
3370       je0 = 0;
3371       do {
3372         WITH = &bridgetable[j - 1];
3373         FORLIM1 = WITH->ie;
3374         for (l = WITH->ib; l <= FORLIM1; l++) {
3375           WITH1 = &chain[l];
3376           for (beta = beta1;
3377                (long)beta <= (long)beta2;
3378                beta = (structure)((long)beta + 1))
3379             P_setcpy(iset[(long)beta - (long)beta1], P_setunion(SET1,
3380                        iset[(long)beta - (long)beta1],
3381                        P_addset(P_expset(SET, 0L),
3382                                 WITH1->ss[(long)beta - (long)symbol])));
3383         }
3384         FORLIM1 = WITH->je;
3385         for (l = WITH->jb; l <= FORLIM1; l++) {
3386           WITH1 = &chain[l];
3387           for (beta = beta1;
3388                (long)beta <= (long)beta2;
3389                beta = (structure)((long)beta + 1))
3390             P_setcpy(jset[(long)beta - (long)beta1], P_setunion(SET1,
3391                        jset[(long)beta - (long)beta1],
3392                        P_addset(P_expset(SET, 0L),
3393                                 WITH1->ss[(long)beta - (long)symbol])));
3394         }
3395         if (WITH->ib < ib0)
3396           ib0 = WITH->ib;
3397         if (WITH->ie > ie0)
3398           ie0 = WITH->ie;
3399         if (WITH->jb < jb0)
3400           jb0 = WITH->jb;
3401         if (WITH->je > je0)
3402           je0 = WITH->je;
3403         j = WITH->towards;
3404       } while (j != 0);
3405       j = i;
3406       if (P_setequal(iset[0], P_addset(P_expset(SET2, 0L), ' ')))
3407         betai = beta1;
3408       else
3409         betai = beta2;
3410       if (P_setequal(jset[0], P_addset(P_expset(SET2, 0L), ' ')))
3411         betaj = beta1;
3412       else
3413         betaj = beta2;
3414       if ((!P_setequal(iset[(long)betai - (long)beta1],
3415                        P_addset(P_expset(SET2, 0L), ' '))) |
3416           (!P_setequal(jset[(long)betaj - (long)beta1],
3417                        P_addset(P_expset(SET3, 0L), ' '))))
3418         if (bVerbose)
3419           printf(" !!! STRAND COLUMN OVERWRITTEN !!!\n");
3420       do {
3421         WITH = &bridgetable[j - 1];
3422         FORLIM1 = WITH->ie;
3423         for (l = WITH->ib; l <= FORLIM1; l++) {
3424           WITH1 = &chain[l];
3425           WITH1->ss[(long)betai - (long)symbol] = WITH->laddername;
3426           if (WITH->btyp == parallel)
3427             WITH1->partner[(long)betai - (long)beta1] = WITH->jb + l - WITH->ib;
3428           else
3429             WITH1->partner[(long)betai - (long)beta1] = WITH->je - l + WITH->ib;
3430         }
3431         FORLIM1 = WITH->je;
3432         for (l = WITH->jb; l <= FORLIM1; l++) {
3433           WITH1 = &chain[l];
3434           WITH1->ss[(long)betaj - (long)symbol] = WITH->laddername;
3435           if (WITH->btyp == parallel)
3436             WITH1->partner[(long)betaj - (long)beta1] = WITH->ib + l - WITH->jb;
3437           else
3438             WITH1->partner[(long)betaj - (long)beta1] = WITH->ie - l + WITH->jb;
3439         }
3440         j = WITH->towards;
3441       } while (j != 0);
3442       if (ib0 == ie0)
3443         cc = 'B';
3444       else
3445         cc = 'E';
3446       for (j = ib0; j <= ie0; j++) {
3447         WITH1 = &chain[j];
3448         if (WITH1->ss[0] != 'E')
3449           WITH1->ss[0] = cc;
3450       }
3451       for (j = jb0; j <= je0; j++) {
3452         WITH1 = &chain[j];
3453         if (WITH1->ss[0] != 'E')
3454           WITH1->ss[0] = cc;
3455       }
3456     }
3457   }
3458   FORLIM = nbridge;
3459   for (j = 0; j < FORLIM; j++) {
3460     WITH = &bridgetable[j];
3461     FORLIM1 = WITH->ie;
3462     for (l = WITH->ib; l <= FORLIM1; l++)
3463       chain[l].sheetlabel = WITH->sheetname;
3464     FORLIM1 = WITH->je;
3465     for (l = WITH->jb; l <= FORLIM1; l++)
3466       chain[l].sheetlabel = WITH->sheetname;
3467   }
3468 }  /* Markstrands */
3469
3470
3471 /***/
3472 /*--------------------------------------------------------------------*/
3473
3474 Static Void Flagbridge()
3475 {
3476   long i, FORLIM;
3477   bridge *WITH;
3478
3479   /***/
3480
3481   for (i = 0; i < MAXBRIDGE; i++) {
3482     WITH = &bridgetable[i];
3483     WITH->ib = 0;
3484     WITH->ie = 0;
3485     WITH->jb = 0;
3486     WITH->je = 0;
3487     WITH->btyp = nobridge;
3488   }
3489   nbridge = 0;
3490   FORLIM = lchain;
3491   for (i = 2; i < FORLIM; i++)
3492     Testbridge(i);
3493   if (nbridge <= 0)
3494     return;
3495   Extendladder();
3496   Sheet();
3497   Markstrands();
3498 }  /* Flagbridge */
3499
3500
3501 /***/
3502
3503 Local Void Flagsymbol()
3504 {
3505   /* FLAGS ALPHA HELICES AND TURNS IN SYMBOL COLUMN */
3506   long i, j, k;
3507   Char cc;
3508   long nhset[9];
3509   structure turn;
3510   boolean empty;
3511   long FORLIM;
3512   backbone *WITH;
3513
3514   P_addset(P_expset(nhset, 0L), '>');
3515   P_addset(nhset, 'X');
3516   FORLIM = lchain - 4;
3517   for (i = 2; i <= FORLIM; i++) {
3518     if (P_inset(chain[i - 1].ss[(long)turn4 - (long)symbol], nhset) &
3519         P_inset(chain[i].ss[(long)turn4 - (long)symbol], nhset)) {
3520       for (j = i; j <= i + 3; j++)
3521         chain[j].ss[0] = 'H';
3522     }
3523   }
3524   FORLIM = lchain - 3;
3525   for (i = 2; i <= FORLIM; i++) {
3526     if (P_inset(chain[i - 1].ss[(long)turn3 - (long)symbol], nhset) &
3527         P_inset(chain[i].ss[(long)turn3 - (long)symbol], nhset)) {
3528       empty = true;
3529       for (j = i; j <= i + 2; j++) {
3530         WITH = &chain[j];
3531         if (WITH->ss[0] != 'G' && WITH->ss[0] != ' ')
3532           empty = false;
3533       }
3534       if (empty) {
3535         for (j = i; j <= i + 2; j++)
3536           chain[j].ss[0] = 'G';
3537       }
3538     }
3539   }
3540   FORLIM = lchain - 5;
3541   for (i = 2; i <= FORLIM; i++) {
3542     if (P_inset(chain[i - 1].ss[(long)turn5 - (long)symbol], nhset) &
3543         P_inset(chain[i].ss[(long)turn5 - (long)symbol], nhset)) {
3544       empty = true;
3545       for (j = i; j <= i + 4; j++) {
3546         WITH = &chain[j];
3547         if (WITH->ss[0] != 'I' && WITH->ss[0] != ' ')
3548           empty = false;
3549       }
3550       if (empty) {
3551         for (j = i; j <= i + 4; j++)
3552           chain[j].ss[0] = 'I';
3553       }
3554     }
3555   }
3556   FORLIM = lchain;
3557   for (i = 2; i < FORLIM; i++) {
3558     WITH = &chain[i];
3559     if (WITH->ss[0] == ' ') {
3560       cc = ' ';
3561       j = 1;
3562       for (turn = turn3;
3563            (long)turn <= (long)turn5;
3564            turn = (structure)((long)turn + 1)) {
3565         j++;
3566         for (k = 1; k <= j; k++) {
3567           if (i > k) {
3568             if (P_inset(chain[i - k].ss[(long)turn - (long)symbol], nhset))
3569               cc = 'T';
3570           }
3571         }
3572       }
3573       if (cc == ' ')
3574         cc = WITH->ss[(long)bend - (long)symbol];
3575       WITH->ss[0] = cc;
3576     }
3577   }
3578 }  /* Flagsymbol */
3579
3580
3581 /***/
3582 /*--------------------------------------------------------------------*/
3583
3584 Static Void Flagturn()
3585 {
3586   long i, j, k;
3587   structure turn;
3588   Char cc;
3589   long FORLIM1;
3590   backbone *WITH;
3591
3592   /***/
3593
3594   k = 2;
3595   cc = '2';
3596   for (turn = turn3; (long)turn <= (long)turn5; turn = (structure)((long)turn + 1)) {
3597     k++;
3598     cc++;
3599     FORLIM1 = lchain - k;
3600     for (i = 1; i <= FORLIM1; i++) {
3601       if (Nochainbreak(i, i + k)) {
3602         if (Testbond(i + k, i)) {
3603           chain[i + k].ss[(long)turn - (long)symbol] = '<';
3604           for (j = 1; j < k; j++) {
3605             WITH = &chain[i + j];
3606             if (WITH->ss[(long)turn - (long)symbol] == ' ')
3607               WITH->ss[(long)turn - (long)symbol] = cc;
3608           }
3609           WITH = &chain[i];
3610           if (WITH->ss[(long)turn - (long)symbol] == '<')
3611             WITH->ss[(long)turn - (long)symbol] = 'X';
3612           else
3613             WITH->ss[(long)turn - (long)symbol] = '>';
3614         }
3615       }
3616     }
3617   }
3618   FORLIM1 = lchain;
3619   for (i = 1; i <= FORLIM1; i++) {
3620     WITH = &chain[i];
3621     if (WITH->kappa != 360.0 && WITH->kappa > 70.0)
3622       WITH->ss[(long)bend - (long)symbol] = 'S';
3623   }
3624   Flagsymbol();
3625 }  /* Flagturn */
3626
3627
3628 /* Local variables for Flagaccess: */
3629 struct LOC_Flagaccess {
3630   long np;
3631   vector p[NFACE];
3632   double wp[NFACE];
3633 } ;
3634
3635 /* Local variables for Polyeder: */
3636 struct LOC_Polyeder {
3637   struct LOC_Flagaccess *LINK;
3638 } ;
3639
3640 /***/
3641
3642 Local Void Triangle(x1, x2, x3, level, LINK)
3643 double *x1, *x2, *x3;
3644 long level;
3645 struct LOC_Polyeder *LINK;
3646 {
3647   long k, level1;
3648   double xnorm;
3649   vector x4, x5, x6;
3650
3651   if (level > 0) {
3652     level1 = level - 1;
3653     for (k = 0; k <= 2; k++) {
3654       x4[k] = x1[k] + x2[k];
3655       x5[k] = x2[k] + x3[k];
3656       x6[k] = x1[k] + x3[k];
3657     }
3658     Norm(x4, &xnorm);
3659     Norm(x5, &xnorm);
3660     Norm(x6, &xnorm);
3661     Triangle(x1, x4, x6, level1, LINK);
3662     Triangle(x4, x2, x5, level1, LINK);
3663     Triangle(x4, x5, x6, level1, LINK);
3664     Triangle(x5, x3, x6, level1, LINK);
3665     return;
3666   }
3667   for (k = 0; k <= 2; k++)
3668     x6[k] = x1[k] + x2[k] + x3[k];
3669   Norm(x6, &xnorm);
3670   LINK->LINK->np++;
3671   VecCopy(LINK->LINK->p[LINK->LINK->np - 1], x6);
3672   Diff(x3, x1, x5);
3673   Diff(x2, x1, x4);
3674   Cross(x5, x4, x6);
3675   Norm(x6, &xnorm);
3676   LINK->LINK->wp[LINK->LINK->np - 1] = xnorm / 2.0;
3677 }  /* Triangle */
3678
3679 /***/
3680
3681 Local Void Polyeder(LINK)
3682 struct LOC_Flagaccess *LINK;
3683 {  /* GENERATES ALL 12 VERTICES OF ICOSAHEDRON */
3684   struct LOC_Polyeder V;
3685   vector v[12];
3686   double a, b;
3687   long i, j, k, level, FORLIM;
3688
3689   /***/
3690
3691   V.LINK = LINK;
3692   k = 0;
3693   a = YVERTEX;
3694   b = ZVERTEX;
3695   for (i = 1; i <= 2; i++) {
3696     a = -a;
3697     for (j = 1; j <= 2; j++) {
3698       b = -b;
3699       k++;
3700       v[k - 1][0] = 0.0;
3701       v[k - 1][1] = a;
3702       v[k - 1][2] = b;
3703       k++;
3704       v[k - 1][0] = b;
3705       v[k - 1][1] = 0.0;
3706       v[k - 1][2] = a;
3707       k++;
3708       v[k - 1][0] = a;
3709       v[k - 1][1] = b;
3710       v[k - 1][2] = 0.0;
3711     }
3712   }
3713   LINK->np = 0;
3714   level = ORDER;
3715   /* GET ALL 20 FACES OF ICOSAHEDRON */
3716   for (i = 0; i <= 9; i++) {   /* FIND INTEGRATION POINTS */
3717     for (j = i + 1; j <= 10; j++) {
3718       if (Distance(v[i], v[j]) < 1.1) {
3719         for (k = j + 1; k <= 11; k++) {
3720           if ((Distance(v[i], v[k]) < 1.1) & (Distance(v[j], v[k]) < 1.1))
3721             Triangle(v[i], v[j], v[k], level, &V);
3722         }
3723       }
3724     }
3725   }
3726   a = 0.0;
3727   FORLIM = LINK->np;
3728   for (i = 0; i < FORLIM; i++)
3729     a += LINK->wp[i];
3730   a = FOURPI / a;
3731   FORLIM = LINK->np;
3732   for (i = 0; i < FORLIM; i++)
3733     LINK->wp[i] *= a;
3734 }  /* Polyeder (enurD idu) */
3735
3736 /* Local variables for Surface: */
3737 struct LOC_Surface {
3738   struct LOC_Flagaccess *LINK;
3739   long nx;
3740   vector x[MAXPACK];
3741   double rx[MAXPACK];
3742 } ;
3743
3744 /***/
3745
3746 Local boolean Step(xx, LINK)
3747 double *xx;
3748 struct LOC_Surface *LINK;
3749 {
3750   long k;
3751   boolean one;
3752   double TEMP;
3753
3754   one = true;
3755   k = 1;
3756   while (k <= LINK->nx && one) {
3757     TEMP = LINK->rx[k - 1] + RWATER;
3758     if (Distsq(xx, LINK->x[k - 1]) < TEMP * TEMP)
3759       one = false;
3760     else
3761       k++;
3762   }
3763   return one;
3764 }  /* Step */
3765
3766 /* Local variables for Liste: */
3767 struct LOC_Liste {
3768   struct LOC_Surface *LINK;
3769 } ;
3770
3771 /***/
3772
3773 Local Void Listentry(xx, yy, d, r, LINK)
3774 double *xx, *yy;
3775 double d, r;
3776 struct LOC_Liste *LINK;
3777 {
3778   vector zz;
3779   double delta;
3780
3781   delta = Distance(xx, yy);
3782   if (delta >= d + r)
3783     return;
3784   if (delta <= EPS)
3785     return;
3786   LINK->LINK->nx++;
3787   if (LINK->LINK->nx > MAXPACK) {
3788     if (bVerbose)
3789       printf(" !!! TABLE OVERFLOW IN FLAGACCESS !!!\n");
3790     longjmp(_JL99, 1);
3791     return;
3792   }
3793   Diff(yy, xx, zz);
3794   VecCopy(LINK->LINK->x[LINK->LINK->nx - 1], zz);
3795   LINK->LINK->rx[LINK->LINK->nx - 1] = r;
3796 }  /* Listentry */
3797
3798 /***/
3799
3800 Local Void Liste(xx, rxx, LINK)
3801 double *xx;
3802 double rxx;
3803 struct LOC_Surface *LINK;
3804 {
3805   struct LOC_Liste V;
3806   long i, k;
3807   double d;
3808   long FORLIM;
3809   backbone *WITH;
3810   long FORLIM1;
3811
3812   /***/
3813
3814   V.LINK = LINK;
3815   LINK->nx = 0;
3816   d = rxx + RWATER + RWATER;
3817   FORLIM = lchain;
3818   for (i = 1; i <= FORLIM; i++) {
3819     if (Nochainbreak(i, i)) {
3820       WITH = &chain[i];
3821       if (Distance(xx, WITH->ca) < d + RESRAD) {
3822         Listentry(xx, WITH->n, d, RN, &V);
3823         Listentry(xx, WITH->ca, d, RCA, &V);
3824         Listentry(xx, WITH->c, d, RC, &V);
3825         Listentry(xx, WITH->o, d, RO, &V);
3826         if (WITH->nsideatoms > 0) {
3827           FORLIM1 = WITH->nsideatoms;
3828           for (k = 0; k < FORLIM1; k++)
3829             Listentry(xx, sidechain[WITH->atompointer + k], d, RSIDEATOM, &V);
3830         }
3831       }
3832     }
3833   }
3834 }  /* Liste */
3835
3836 /***/
3837
3838 Local double Surface(xatom, ratom, LINK)
3839 double *xatom;
3840 double ratom;
3841 struct LOC_Flagaccess *LINK;
3842 {
3843   struct LOC_Surface V;
3844   long i, j;
3845   double f, radius;
3846   vector xx;
3847   long FORLIM;
3848
3849   /***/
3850
3851   V.LINK = LINK;
3852   Liste(xatom, ratom, &V);
3853   radius = ratom + RWATER;
3854   f = 0.0;
3855   FORLIM = LINK->np;
3856   for (i = 0; i < FORLIM; i++) {
3857     for (j = 0; j <= 2; j++)
3858       xx[j] = LINK->p[i][j] * radius;
3859     if (Step(xx, &V))
3860       f += LINK->wp[i];
3861   }
3862   return (radius * radius * f);
3863 }  /* Surface */
3864
3865
3866 /***/
3867 /*--------------------------------------------------------------------*/
3868
3869 Static Void Flagaccess()
3870 {
3871   struct LOC_Flagaccess V;
3872   long i, k;
3873   double f;
3874   long FORLIM;
3875   backbone *WITH;
3876   long FORLIM1;
3877
3878   /***/
3879
3880   Polyeder(&V);
3881   FORLIM = lchain;
3882   for (i = 1; i <= FORLIM; i++) {
3883     if (Nochainbreak(i, i)) {
3884       WITH = &chain[i];
3885       f = Surface(WITH->n, RN, &V) + Surface(WITH->ca, RCA, &V) +
3886           Surface(WITH->c, RC, &V) + Surface(WITH->o, RO, &V);
3887       if (WITH->nsideatoms > 0) {
3888         FORLIM1 = WITH->nsideatoms;
3889         for (k = 0; k < FORLIM1; k++)
3890           f += Surface(sidechain[WITH->atompointer + k], RSIDEATOM, &V);
3891       }
3892       WITH->access = (long)floor(f + 0.5);
3893     }
3894   }
3895 }  /* Flagaccess */
3896
3897
3898 /***/
3899
3900 Local Void Statistics()
3901 {
3902   long i, j, k, nchain, nres, nhbond, lhelix;
3903   bridgetyp b;
3904   Char cc;
3905   double Surface;
3906   long nhbturn[11];
3907   long ladderset[MAXBRIDGE / 32 + 2];
3908   long hbridge[(long)antiparallel - (long)parallel + 1];
3909   long helixhist[MAXHIST], sheethist[MAXHIST];
3910   long betahist[(long)antiparallel - (long)parallel + 1][MAXHIST];
3911   long FORLIM, FORLIM1;
3912   backbone *WITH;
3913   bridge *WITH1;
3914   long SET[11];
3915   long SET1[257];
3916
3917   lhelix = 0;
3918   nhbond = 0;
3919   nchain = 0;
3920   nres = 0;
3921   for (i = 0; i < MAXHIST; i++) {
3922     for (b = parallel; (long)b <= (long)antiparallel; b = (bridgetyp)((long)b + 1))
3923       betahist[(long)b - (long)parallel][i] = 0;
3924     helixhist[i] = 0;
3925     sheethist[i] = 0;
3926   }
3927   Surface = 0.0;
3928   for (k = 0; k <= 10; k++)
3929     nhbturn[k] = 0;
3930   for (b = parallel; (long)b <= (long)antiparallel; b = (bridgetyp)((long)b + 1))
3931     hbridge[(long)b - (long)parallel] = 0;
3932   FORLIM = lchain;
3933   for (i = 0; i <= FORLIM; i++) {
3934     WITH = &chain[i];
3935     if (Nochainbreak(i, i)) {
3936       nres++;
3937       Surface += WITH->access;
3938       for (j = 0; j <= 1; j++) {
3939         if (WITH->donor[j].energy < HBHIGH) {
3940           nhbond++;
3941           k = WITH->donor[j].residue - i;
3942           if (labs(k) < 6)
3943             nhbturn[k + 5]++;
3944         }
3945       }
3946     } else
3947       nchain++;
3948     if (WITH->ss[0] == 'H')
3949       lhelix++;
3950     else if (lhelix > 0) {
3951       if (lhelix > MAXHIST)
3952         lhelix = MAXHIST;
3953       helixhist[lhelix - 1]++;
3954       lhelix = 0;
3955     }
3956   }
3957   if (nbridge > 0) {
3958     FORLIM = nbridge;
3959     for (i = 1; i <= FORLIM; i++) {
3960       WITH1 = &bridgetable[i - 1];
3961       hbridge[(long)WITH1->btyp - (long)parallel] += WITH1->ie - WITH1->ib + 2;
3962       if (WITH1->from == 0) {
3963         j = i;
3964         k = 0;
3965         do {
3966           k += bridgetable[j - 1].ie - bridgetable[j - 1].ib + 1;
3967           j = bridgetable[j - 1].towards;
3968         } while (j != 0);
3969         if (k > MAXHIST)
3970           k = MAXHIST;
3971         betahist[(long)WITH1->btyp - (long)parallel][k - 1]++;
3972       }
3973     }
3974   }
3975   if (nbridge > 0) {
3976     P_expset(ladderset, 0L);
3977     FORLIM = nbridge;
3978     for (i = 1; i <= FORLIM; i++)
3979       P_addset(ladderset, (int)i);
3980     FORLIM = nbridge;
3981     for (i = 1; i <= FORLIM; i++) {
3982       WITH1 = &bridgetable[i - 1];
3983       if ((WITH1->from == 0) & P_inset((int)i, ladderset)) {
3984         if (!P_setequal(P_addset(P_expset(SET1, 0L), (int)i), WITH1->linkset) ||
3985             WITH1->ie > WITH1->ib) {
3986           k = 0;
3987           FORLIM1 = nbridge;
3988           for (j = 1; j <= FORLIM1; j++) {
3989             if ((bridgetable[j - 1].from == 0) & P_inset((int)j, WITH1->linkset))
3990               k++;
3991           }
3992           sheethist[k - 1]++;
3993         }
3994         P_setdiff(ladderset, ladderset, WITH1->linkset);
3995       }
3996     }
3997   }
3998   if (nres == 0)
3999     return;
4000   fprintf(tapeout,
4001     "%5ld%3ld%3ld%3ld%3ld TOTAL NUMBER OF RESIDUES, NUMBER OF CHAINS, NUMBER OF SS-BRIDGES(TOTAL,INTRACHAIN,INTERCHAIN)                .\n",
4002     nres, nchain, nssinter + nssintra, nssintra, nssinter);
4003   fprintf(tapeout,
4004     "%8.1f   ACCESSIBLE SURFACE OF PROTEIN (ANGSTROM**2)                                                                         .\n",
4005     Surface);
4006   fprintf(tapeout,
4007     "%5ld%5.1f   TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(J)  , SAME NUMBER PER 100 RESIDUES                              .\n",
4008     nhbond, 100.0 * nhbond / nres);
4009   i = hbridge[0];
4010   j = hbridge[(long)antiparallel - (long)parallel];
4011   fprintf(tapeout,
4012     "%5ld%5.1f   TOTAL NUMBER OF HYDROGEN BONDS IN     PARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES                              .\n",
4013     i, 100.0 * i / nres);
4014   fprintf(tapeout,
4015     "%5ld%5.1f   TOTAL NUMBER OF HYDROGEN BONDS IN ANTIPARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES                              .\n",
4016     j, 100.0 * j / nres);
4017   for (i = -5; i <= 5; i++) {
4018     if (i < 0)
4019       cc = '-';
4020     else
4021       cc = '+';
4022     k = labs(i);
4023     fprintf(tapeout,
4024       "%5ld%5.1f   TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(I%c%ld), SAME NUMBER PER 100 RESIDUES                              .\n",
4025       nhbturn[i + 5], 100.0 * nhbturn[i + 5] / nres, cc, k);
4026   }
4027   for (i = 1; i <= MAXHIST; i++)
4028     fprintf(tapeout, "%3ld", i);
4029   fprintf(tapeout, "     *** HISTOGRAMS OF ***           .\n");
4030   for (i = 0; i < MAXHIST; i++)
4031     fprintf(tapeout, "%3ld", helixhist[i]);
4032   fprintf(tapeout, "    RESIDUES PER ALPHA HELIX         .\n");
4033   for (i = 0; i < MAXHIST; i++)
4034     fprintf(tapeout, "%3ld", betahist[0][i]);
4035   fprintf(tapeout, "    PARALLEL BRIDGES PER LADDER      .\n");
4036   for (i = 0; i < MAXHIST; i++)
4037     fprintf(tapeout, "%3ld", betahist[(long)antiparallel - (long)parallel][i]);
4038   fprintf(tapeout, "    ANTIPARALLEL BRIDGES PER LADDER  .\n");
4039   for (i = 0; i < MAXHIST; i++)
4040     fprintf(tapeout, "%3ld", sheethist[i]);
4041   fprintf(tapeout, "    LADDERS PER SHEET                .\n");
4042 }  /* Statistics */
4043
4044 /***/
4045
4046 Local Void Writehb(i, hb)
4047 long i;
4048 hydrogenbond hb;
4049 {
4050   double e;
4051
4052   if (hb.residue != 0)
4053     hb.residue -= i;
4054   e = hb.energy / 1000.0;
4055   fprintf(tapeout, "%4ld,%4.1f", hb.residue, e);
4056 }  /* Writehb */
4057
4058
4059 /***/
4060 /*--------------------------------------------------------------------*/
4061
4062 Static Void Printout()
4063 {
4064   long i, j;
4065   structure s;
4066   double phi, psi, tco;
4067   long FORLIM;
4068   backbone *WITH;
4069
4070   /***/
4071
4072   Statistics();
4073   fprintf(tapeout,
4074     "  #  RESIDUE AA STRUCTURE BP1 BP2  ACC   N-H-->O  O-->H-N  N-H-->O  O-->H-N    TCO  KAPPA ALPHA  PHI   PSI    X-CA   Y-CA   Z-CA \n");
4075   FORLIM = lchain;
4076   for (i = 1; i <= FORLIM; i++) {
4077     WITH = &chain[i];
4078     fprintf(tapeout, "%5ld ", i);
4079     for (j = 0; j <= 5; j++)
4080       putc(WITH->aaident[j], tapeout);
4081     fprintf(tapeout, " %c  %c ", WITH->aa, WITH->ss[0]);
4082     for (s = turn3; (long)s <= (long)beta2; s = (structure)((long)s + 1))
4083       putc(WITH->ss[(long)s - (long)symbol], tapeout);
4084     for (s = beta1; (long)s <= (long)beta2; s = (structure)((long)s + 1))
4085       fprintf(tapeout, "%4ld", WITH->partner[(long)s - (long)beta1]);
4086     fprintf(tapeout, "%c%4ld ", WITH->sheetlabel, WITH->access);
4087     for (j = 0; j <= 1; j++) {
4088       Writehb(i, WITH->acceptor[j]);
4089       Writehb(i, WITH->donor[j]);
4090     }
4091     phi = 360.0;
4092     psi = 360.0;
4093     tco = 0.0;
4094     if (Nochainbreak(i - 1, i)) {
4095       phi = Dihedralangle(chain[i - 1].c, WITH->n, WITH->ca, WITH->c);
4096       tco = Cosangle(WITH->c, WITH->o, chain[i - 1].c, chain[i - 1].o);
4097     }
4098     if (Nochainbreak(i, i + 1))
4099       psi = Dihedralangle(WITH->n, WITH->ca, WITH->c, chain[i + 1].n);
4100     fprintf(tapeout, "%8.3f%6.1f%6.1f%6.1f%6.1f%7.1f%7.1f%7.1f\n",
4101             tco, WITH->kappa, WITH->alpha, phi, psi, WITH->ca[0], WITH->ca[1],
4102             WITH->ca[2]);
4103   }
4104 }  /* Printout */
4105
4106 Static Void Usage()
4107 {
4108   fprintf(stderr,"Usage: dssp [-na] pdb_file dssp_file\n");
4109   fprintf(stderr,"the -na flag disables the calculation of accessible surface\n");
4110 }
4111
4112 /***/
4113 /*--------------------------------------------------------------------*/
4114
4115 void printit()
4116 {
4117   printf(" \n");
4118   printf("                           DSSP\n");
4119   printf("            by Wolfgang Kabsch and Chris Sander\n");
4120   printf("\n");
4121   printf("Defines secondary structure and solvent exposure of proteins from\n");
4122   printf("atomic coordinates as given in Protein Data Bank format.   \n");
4123   printf("\n");
4124   printf("_________________________________________________________________________\n");
4125   printf("This version licensed to Ethan Benatan at Univ_Pittsburgh               \n");
4126   printf("for academic purposes.                                                 \n");
4127   printf("Do not redistribute. \n");
4128   printf("\n");
4129   printf("Commercial licenses available on request.                                \n");
4130   printf("\n");
4131   printf("Copyright by Wolfgang Kabsch and Chris Sander, 1983, 1985, 1988. \n");
4132   printf("Fax: +49-6221-387 306\n");
4133   printf("\n");
4134   printf("Algorithm version October 1988. Refer to Biopolymers 22(1983) 2577-2637.\n");
4135   printf("Do report errors if you find any.\n");
4136   printf("\n");
4137   printf("Email: Sander@embl-heidelberg.de \n");
4138   printf("       Kabsch@embl-heidelberg.de \n");
4139   printf("\n");
4140   printf("_________________________________________________________________________\n");
4141   printf("Related databases and datasets available from the Protein\n");
4142   printf("Design Group at EMBL via anonymous ftp from embl-heidelberg.de:\n");
4143   printf("\n");
4144   printf("pdb_select   Representative set of protein structures.\n");
4145   printf("             By U. Hobohm, C. Sander, M. Scharf and R. Schneider.\n");
4146   printf("             See Protein Science 1, 409-417.\n");
4147   printf("DSSP         Dictionary of secondary structures of proteins. \n");
4148   printf("HSSP         Database of sequence-homology derived protein families.\n");
4149   printf("             By C. Sander and R.Schneider.\n");
4150   printf("             See Proteins 9, 56-68 (1991).\n");
4151   printf("FSSP         Database of protein structure families with \n");
4152   printf("             common folding motifs. \n");
4153   printf("             L.Holm, C. Ouzounis, C. Sander, G.Tuparev, G. Vriend\n");
4154   printf("             See Protein Science, in the press (1992).\n");
4155   printf("In the XSSP databases, there is one dataset for each unique or\n");
4156   printf("             representative PDB protein, e.g., 1PPT.HSSP etc.\n");
4157   printf("\n");
4158   printf("Restrictions:Commercial users must apply for a license. \n");
4159   printf("             Not to be used for classified research.\n");
4160 }
4161
4162 void dssp_main(int bDoAcc, int bSetVerbose)
4163 {
4164   int tt; /*TIMELOCK*/
4165   PASCAL_MAIN(0,NULL);
4166   if (setjmp(_JL99))
4167     goto _L99;
4168
4169   bVerbose=bSetVerbose;
4170
4171   tt=time(0); /*TIMELOCK*/
4172
4173   lchain = 0;
4174   Inputcoordinates(&lchain);
4175   if (!Nochainbreak(1L, lchain))
4176     printf(" !!! POLYPEPTIDE CHAIN INTERRUPTED !!!\n");
4177 /*   printf("INPUTCOORDINATES DONE%12ld\n", lchain); */
4178   Flagssbonds();
4179 /*   printf("FLAGSSBONDS DONE\n"); */
4180   Flagchirality();
4181 /*   printf("FLAGCHIRALITY DONE\n"); */
4182   Flaghydrogenbonds();
4183 /*   printf("FLAGHYDROGENBONDS DONE\n"); */
4184   Flagbridge();
4185 /*   printf("FLAGBRIDGE DONE\n"); */
4186   Flagturn();
4187 /*   printf("FLAGTURN DONE\n"); */
4188   if (bDoAcc) {
4189     Flagaccess();
4190 /*     printf("FLAGACCESS DONE\n"); */
4191   }
4192   Printout();
4193 /*   printf("PRINTOUT DONE\n"); */
4194 _L99:
4195   /*fclose(tapein);
4196     fclose(tapeout);*/
4197   ;
4198 }  /* END OF PROGRAM DSSP */
4199