src/local/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 #ifndef HAS_STDLIB
 606 extern Anyptr   malloc      PP( (size_t) );
 607 extern Void     free        PP( (Anyptr) );
 608 #endif
 609
 610 extern int      _OutMem     PV();
 611 extern int      _CaseCheck  PV();
 612 extern int      _NilCheck   PV();
 613 extern int      _Escape     PP( (int) );
 614 extern int      _EscIO      PP( (int) );
 615
 616 extern long     ipow        PP( (long, long) );
 617 extern Char    *strsub      PP( (Char *, Char *, int, int) );
 618 extern Char    *strltrim    PP( (Char *) );
 619 extern Char    *strrtrim    PP( (Char *) );
 620 extern Char    *strrpt      PP( (Char *, Char *, int) );
 621 extern Char    *strpad      PP( (Char *, Char *, int, int) );
 622 extern int      strpos2     PP( (Char *, Char *, int) );
 623 extern long     memavail    PV();
 624 extern int      P_peek      PP( (FILE *) );
 625 extern int      P_eof       PP( (FILE *) );
 626 extern int      P_eoln      PP( (FILE *) );
 627 extern Void     P_readpaoc  PP( (FILE *, Char *, int) );
 628 extern Void     P_readlnpaoc PP( (FILE *, Char *, int) );
 629 extern long     P_maxpos    PP( (FILE *) );
 630 extern Char    *P_trimname  PP( (Char *, int) );
 631 extern long    *P_setunion  PP( (long *, long *, long *) );
 632 extern long    *P_setint    PP( (long *, long *, long *) );
 633 extern long    *P_setdiff   PP( (long *, long *, long *) );
 634 extern long    *P_setxor    PP( (long *, long *, long *) );
 635 extern int      P_inset     PP( (unsigned, long *) );
 636 extern int      P_setequal  PP( (long *, long *) );
 637 extern int      P_subset    PP( (long *, long *) );
 638 extern long    *P_addset    PP( (long *, unsigned) );
 639 extern long    *P_addsetr   PP( (long *, unsigned, unsigned) );
 640 extern long    *P_remset    PP( (long *, unsigned) );
 641 extern long    *P_setcpy    PP( (long *, long *) );
 642 extern long    *P_expset    PP( (long *, long) );
 643 extern long     P_packset   PP( (long *) );
 644 extern int      P_getcmdline PP( (int l, int h, Char *line) );
 645 extern Void     TimeStamp   PP( (int *Day, int *Month, int *Year,
 646                                  int *Hour, int *Min, int *Sec) );
 647 extern Void     P_sun_argv  PP( (char *, int, int) );
 648
 649
 650 /* I/O error handling */
 651 #define _CHKIO(cond,ior,val,def)  ((cond) ? P_ioresult=0,(val)  \
 652                                           : P_ioresult=(ior),(def))
 653 #define _SETIO(cond,ior)          (P_ioresult = (cond) ? 0 : (ior))
 654
 655 /* Following defines are suitable for the HP Pascal operating system */
 656 #define FileNotFound     10
 657 #define FileNotOpen      13
 658 #define FileWriteError   38
 659 #define BadInputFormat   14
 660 #define EndOfFile        30
 661
 662 /* Creating temporary files */
 663 #if (defined(BSD) || defined(NO_TMPFILE)) && !defined(HAVE_TMPFILE)
 664 # define tmpfile()  (fopen(tmpnam(NULL), "w+"))
 665 #endif
 666
 667 /* File buffers */
 668 #define FILEBUF(f,sc,type) sc int __CAT__(f,_BFLAGS);   \
 669                            sc type __CAT__(f,_BUFFER)
 670 #define FILEBUFNC(f,type)  int __CAT__(f,_BFLAGS);   \
 671                            type __CAT__(f,_BUFFER)
 672
 673 #define RESETBUF(f,type)   (__CAT__(f,_BFLAGS) = 1)
 674 #define SETUPBUF(f,type)   (__CAT__(f,_BFLAGS) = 0)
 675
 676 #define GETFBUF(f,type)    (*((__CAT__(f,_BFLAGS) == 1 &&   \
 677                                ((__CAT__(f,_BFLAGS) = 2),   \
 678                                 fread(&__CAT__(f,_BUFFER),  \
 679                                       sizeof(type),1,(f)))),\
 680                               &__CAT__(f,_BUFFER)))
 681 #define AGETFBUF(f,type)   ((__CAT__(f,_BFLAGS) == 1 &&   \
 682                              ((__CAT__(f,_BFLAGS) = 2),   \
 683                               fread(__CAT__(f,_BUFFER),  \
 684                                     sizeof(type),1,(f)))),\
 685                             __CAT__(f,_BUFFER))
 686
 687 #define PUTFBUF(f,type,v)  (GETFBUF(f,type) = (v))
 688 #define CPUTFBUF(f,v)      (PUTFBUF(f,char,v))
 689 #define APUTFBUF(f,type,v) (memcpy(AGETFBUF(f,type), (v),  \
 690                                    sizeof(__CAT__(f,_BUFFER))))
 691
 692 #define GET(f,type)        (__CAT__(f,_BFLAGS) == 1 ?   \
 693                             fread(&__CAT__(f,_BUFFER),sizeof(type),1,(f)) :  \
 694                             (__CAT__(f,_BFLAGS) = 1))
 695
 696 #define PUT(f,type)        (fwrite(&__CAT__(f,_BUFFER),sizeof(type),1,(f)),  \
 697                             (__CAT__(f,_BFLAGS) = 0))
 698 #define CPUT(f)            (PUT(f,char))
 699
 700 #define BUFEOF(f)          (__CAT__(f,_BFLAGS) != 2 && P_eof(f))
 701 #define BUFFPOS(f)         (ftell(f) - (__CAT__(f,_BFLAGS) == 2))
 702
 703 typedef struct {
 704     FILE *f;
 705     FILEBUFNC(f,Char);
 706     Char name[_FNSIZE];
 707 } _TEXT;
 708
 709 /* Memory allocation */
 710 #ifdef __GCC__
 711 # define Malloc(n)  (malloc(n) ?: (Anyptr)_OutMem())
 712 #else
 713 extern Anyptr __MallocTemp__;
 714 # define Malloc(n)  ((__MallocTemp__ = malloc(n)) ? __MallocTemp__ : (Anyptr)_OutMem())
 715 #endif
 716 #define FreeR(p)    (free((Anyptr)(p)))    /* used if arg is an rvalue */
 717 #define Free(p)     (free((Anyptr)(p)), (p)=NULL)
 718
 719 /* sign extension */
 720 #define SEXT(x,n)   ((x) | -(((x) & (1L<<((n)-1))) << 1))
 721
 722 /* packed arrays */   /* BEWARE: these are untested! */
 723 #define P_getbits_UB(a,i,n,L)   ((int)((a)[(i)>>(L)-(n)] >>   \
 724                                        (((~(i))&((1<<(L)-(n))-1)) << (n)) &  \
 725                                        (1<<(1<<(n)))-1))
 726
 727 #define P_getbits_SB(a,i,n,L)   ((int)((a)[(i)>>(L)-(n)] <<   \
 728                                        (16 - ((((~(i))&((1<<(L)-(n))-1))+1) <<\
 729                                               (n)) >> (16-(1<<(n))))))
 730
 731 #define P_putbits_UB(a,i,x,n,L) ((a)[(i)>>(L)-(n)] |=   \
 732                                  (x) << (((~(i))&((1<<(L)-(n))-1)) << (n)))
 733
 734 #define P_putbits_SB(a,i,x,n,L) ((a)[(i)>>(L)-(n)] |=   \
 735                                  ((x) & (1<<(1<<(n)))-1) <<   \
 736                                  (((~(i))&((1<<(L)-(n))-1)) << (n)))
 737
 738 #define P_clrbits_B(a,i,n,L)    ((a)[(i)>>(L)-(n)] &=   \
 739                                  ~( ((1<<(1<<(n)))-1) <<   \
 740                                    (((~(i))&((1<<(L)-(n))-1)) << (n))) )
 741
 742 /* small packed arrays */
 743 #define P_getbits_US(v,i,n)     ((int)((v) >> ((i)<<(n)) & (1<<(1<<(n)))-1))
 744 #define P_getbits_SS(v,i,n)     ((int)((long)(v) << (SETBITS - (((i)+1) << (n))) >> (SETBITS-(1<<(n)))))
 745 #define P_putbits_US(v,i,x,n)   ((v) |= (x) << ((i) << (n)))
 746 #define P_putbits_SS(v,i,x,n)   ((v) |= ((x) & (1<<(1<<(n)))-1) << ((i)<<(n)))
 747 #define P_clrbits_S(v,i,n)      ((v) &= ~( ((1<<(1<<(n)))-1) << ((i)<<(n)) ))
 748
 749 #define P_max(a,b)   ((a) > (b) ? (a) : (b))
 750 #define P_min(a,b)   ((a) < (b) ? (a) : (b))
 751
 752
 753 /* Fix ANSI-isms */
 754
 755 #ifdef LACK_LABS
 756 # ifndef labs
 757 #  define labs  my_labs
 758    extern long my_labs PP( (long) );
 759 # endif
 760 #endif
 761
 762 #ifdef LACK_MEMMOVE
 763 # ifndef memmove
 764 #  define memmove  my_memmove
 765    extern Anyptr my_memmove PP( (Anyptr, Const Anyptr, size_t) );
 766 # endif
 767 #endif
 768
 769 #ifdef LACK_MEMCPY
 770 # ifndef memcpy
 771 #  define memcpy  my_memcpy
 772    extern Anyptr my_memcpy PP( (Anyptr, Const Anyptr, size_t) );
 773 # endif
 774 # ifndef memcmp
 775 #  define memcmp  my_memcmp
 776    extern int my_memcmp PP( (Const Anyptr, Const Anyptr, size_t) );
 777 # endif
 778 # ifndef memset
 779 #  define memset  my_memset
 780    extern Anyptr my_memset PP( (Anyptr, int, size_t) );
 781 # endif
 782 #endif
 783
 784 /* Fix toupper/tolower on Suns and other stupid BSD systems */
 785 #ifdef toupper
 786 # undef toupper
 787 # undef tolower
 788 # define toupper(c)   my_toupper(c)
 789 # define tolower(c)   my_tolower(c)
 790 #endif
 791
 792 #ifndef _toupper
 793 # if 'A' == 65 && 'a' == 97
 794 #  define _toupper(c)  ((c)-'a'+'A')
 795 #  define _tolower(c)  ((c)-'A'+'a')
 796 # else
 797 #  ifdef toupper
 798 #   undef toupper   /* hope these are shadowing real functions, */
 799 #   undef tolower   /* because my_toupper calls _toupper! */
 800 #  endif
 801 #  define _toupper(c)  toupper(c)
 802 #  define _tolower(c)  tolower(c)
 803 # endif
 804 #endif
 805
 806
 807 #endif    /* P2C_H */
 808
 809
 810
 811 /* End. */
 812
 813
 814
 815
 816 /************************************************************
 817 *                                                           *
 818 *          p2clib.c                                           *
 819 *                                                           *
 820 ************************************************************/
 821
 822 /* Run-time library for use with "p2c", the Pascal to C translator */
 823
 824 /* "p2c"  Copyright (C) 1989, 1990, 1991 Free Software Foundation.
 825  * By Dave Gillespie, daveg@csvax.cs.caltech.edu.  Version --VERSION--.
 826  * This file may be copied, modified, etc. in any way.  It is not restricted
 827  * by the licence agreement accompanying p2c itself.
 828  */
 829
 830
 831
 832  /* #include "p2c.h" */
 833
 834
 835 #ifndef NO_TIME
 836 # include <time.h>
 837 #endif
 838
 839
 840 #define Isspace(c)  isspace(c)      /* or "((c) == ' ')" if preferred */
 841
 842
 843
 844
 845 int P_argc;
 846 char **P_argv;
 847
 848 short P_escapecode;
 849 int P_ioresult;
 850
 851 long EXCP_LINE;    /* Used by Pascal workstation system */
 852
 853 Anyptr __MallocTemp__;
 854
 855 __p2c_jmp_buf *__top_jb;
 856
 857
 858
 859
 860 void PASCAL_MAIN(argc, argv)
 861 int argc;
 862 char **argv;
 863 {
 864     P_argc = argc;
 865     P_argv = argv;
 866     __top_jb = NULL;
 867
 868 #ifdef LOCAL_INIT
 869     LOCAL_INIT();
 870 #endif
 871 }
 872
 873
 874
 875
 876
 877 /* In case your system lacks these... */
 878
 879 long my_labs(x)
 880 long x;
 881 {
 882     return((x > 0) ? x : -x);
 883 }
 884
 885
 886 #ifdef __STDC__
 887 Anyptr my_memmove(Anyptr d, Const Anyptr s, size_t n)
 888 #else
 889 Anyptr my_memmove(d, s, n)
 890 Anyptr d, s;
 891 register long n;
 892 #endif
 893 {
 894     register char *dd = (char *)d, *ss = (char *)s;
 895     if (dd < ss || dd - ss >= n) {
 896         memcpy(dd, ss, n);
 897     } else if (n > 0) {
 898         dd += n;
 899         ss += n;
 900         while (--n >= 0)
 901             *--dd = *--ss;
 902     }
 903     return d;
 904 }
 905
 906
 907 #ifdef __STDC__
 908 Anyptr my_memcpy(Anyptr d, Const Anyptr s, size_t n)
 909 #else
 910 Anyptr my_memcpy(d, s, n)
 911 Anyptr d, s;
 912 register long n;
 913 #endif
 914 {
 915     register char *ss = (char *)s, *dd = (char *)d;
 916     while (--n >= 0)
 917         *dd++ = *ss++;
 918     return d;
 919 }
 920
 921 #ifdef __STDC__
 922 int my_memcmp(Const Anyptr s1, Const Anyptr s2, size_t n)
 923 #else
 924 int my_memcmp(s1, s2, n)
 925 Anyptr s1, s2;
 926 register long n;
 927 #endif
 928 {
 929     register char *a = (char *)s1, *b = (char *)s2;
 930     register int i;
 931     while (--n >= 0)
 932         if ((i = (*a++) - (*b++)) != 0)
 933             return i;
 934     return 0;
 935 }
 936
 937 #ifdef __STDC__
 938 Anyptr my_memset(Anyptr d, int c, size_t n)
 939 #else
 940 Anyptr my_memset(d, c, n)
 941 Anyptr d;
 942 register int c;
 943 register long n;
 944 #endif
 945 {
 946     register char *dd = (char *)d;
 947     while (--n >= 0)
 948         *dd++ = c;
 949     return d;
 950 }
 951
 952
 953 int my_toupper(c)
 954 int c;
 955 {
 956     if (islower(c))
 957         return _toupper(c);
 958     else
 959         return c;
 960 }
 961
 962
 963 int my_tolower(c)
 964 int c;
 965 {
 966     if (isupper(c))
 967         return _tolower(c);
 968     else
 969         return c;
 970 }
 971
 972
 973
 974
 975 long ipow(a, b)
 976 long a, b;
 977 {
 978     long v;
 979
 980     if (a == 0 || a == 1)
 981         return a;
 982     if (a == -1)
 983         return (b & 1) ? -1 : 1;
 984     if (b < 0)
 985         return 0;
 986     if (a == 2)
 987         return 1 << b;
 988     v = (b & 1) ? a : 1;
 989     while ((b >>= 1) > 0) {
 990         a *= a;
 991         if (b & 1)
 992             v *= a;
 993     }
 994     return v;
 995 }
 996
 997
 998
 999
1000 /* Common string functions: */
1001
1002 /* Store in "ret" the substring of length "len" starting from "pos" (1-based).
1003    Store a shorter or null string if out-of-range.  Return "ret". */
1004
1005 char *strsub(ret, s, pos, len)
1006 register char *ret, *s;
1007 register int pos, len;
1008 {
1009     register char *s2;
1010
1011     if (--pos < 0 || len <= 0) {
1012         *ret = 0;
1013         return ret;
1014     }
1015     while (pos > 0) {
1016         if (!*s++) {
1017             *ret = 0;
1018             return ret;
1019         }
1020         pos--;
1021     }
1022     s2 = ret;
1023     while (--len >= 0) {
1024         if (!(*s2++ = *s++))
1025             return ret;
1026     }
1027     *s2 = 0;
1028     return ret;
1029 }
1030
1031
1032 /* Return the index of the first occurrence of "pat" as a substring of "s",
1033    starting at index "pos" (1-based).  Result is 1-based, 0 if not found. */
1034
1035 int strpos2(s, pat, pos)
1036 char *s;
1037 register char *pat;
1038 register int pos;
1039 {
1040     register char *cp, ch;
1041     register int slen;
1042
1043     if (--pos < 0)
1044         return 0;
1045     slen = strlen(s) - pos;
1046     cp = s + pos;
1047     if (!(ch = *pat++))
1048         return 0;
1049     pos = strlen(pat);
1050     slen -= pos;
1051     while (--slen >= 0) {
1052         if (*cp++ == ch && !strncmp(cp, pat, pos))
1053             return cp - s;
1054     }
1055     return 0;
1056 }
1057
1058
1059 /* Case-insensitive version of strcmp. */
1060
1061 int strcicmp(s1, s2)
1062 register char *s1, *s2;
1063 {
1064     register unsigned char c1, c2;
1065
1066     while (*s1) {
1067         if (*s1++ != *s2++) {
1068             if (!s2[-1])
1069                 return 1;
1070             c1 = toupper(s1[-1]);
1071             c2 = toupper(s2[-1]);
1072             if (c1 != c2)
1073                 return c1 - c2;
1074         }
1075     }
1076     if (*s2)
1077         return -1;
1078     return 0;
1079 }
1080
1081
1082
1083
1084 /* HP and Turbo Pascal string functions: */
1085
1086 /* Trim blanks at left end of string. */
1087
1088 char *strltrim(s)
1089 register char *s;
1090 {
1091     while (Isspace(*s++)) ;
1092     return s - 1;
1093 }
1094
1095
1096 /* Trim blanks at right end of string. */
1097
1098 char *strrtrim(s)
1099 register char *s;
1100 {
1101     register char *s2 = s;
1102
1103     if (!*s)
1104         return s;
1105     while (*++s2) ;
1106     while (s2 > s && Isspace(*--s2))
1107         *s2 = 0;
1108     return s;
1109 }
1110
1111
1112 /* Store in "ret" "num" copies of string "s".  Return "ret". */
1113
1114 char *strrpt(ret, s, num)
1115 char *ret;
1116 register char *s;
1117 register int num;
1118 {
1119     register char *s2 = ret;
1120     register char *s1;
1121
1122     while (--num >= 0) {
1123         s1 = s;
1124         while ((*s2++ = *s1++)) ;
1125         s2--;
1126     }
1127     return ret;
1128 }
1129
1130
1131 /* Store in "ret" string "s" with enough pad chars added to reach "size". */
1132
1133 char *strpad(ret, s, padchar, num)
1134 char *ret;
1135 register char *s;
1136 register int padchar, num;
1137 {
1138     register char *d = ret;
1139
1140     if (s == d) {
1141         while (*d++) ;
1142     } else {
1143         while ((*d++ = *s++)) ;
1144     }
1145     num -= (--d - ret);
1146     while (--num >= 0)
1147         *d++ = padchar;
1148     *d = 0;
1149     return ret;
1150 }
1151
1152
1153 /* Copy the substring of length "len" from index "spos" of "s" (1-based)
1154    to index "dpos" of "d", lengthening "d" if necessary.  Length and
1155    indices must be in-range. */
1156
1157 void strmove(len, s, spos, d, dpos)
1158 register char *s, *d;
1159 register int len, spos, dpos;
1160 {
1161     s += spos - 1;
1162     d += dpos - 1;
1163     while (*d && --len >= 0)
1164         *d++ = *s++;
1165     if (len > 0) {
1166         while (--len >= 0)
1167             *d++ = *s++;
1168         *d = 0;
1169     }
1170 }
1171
1172
1173 /* Delete the substring of length "len" at index "pos" from "s".
1174    Delete less if out-of-range. */
1175
1176 void strdelete(s, pos, len)
1177 register char *s;
1178 register int pos, len;
1179 {
1180     register int slen;
1181
1182     if (--pos < 0)
1183         return;
1184     slen = strlen(s) - pos;
1185     if (slen <= 0)
1186         return;
1187     s += pos;
1188     if (slen <= len) {
1189         *s = 0;
1190         return;
1191     }
1192     while ((*s = s[len])) s++;
1193 }
1194
1195
1196 /* Insert string "src" at index "pos" of "dst". */
1197
1198 void strinsert(src, dst, pos)
1199 register char *src, *dst;
1200 register int pos;
1201 {
1202     register int slen, dlen;
1203
1204     if (--pos < 0)
1205         return;
1206     dlen = strlen(dst);
1207     dst += dlen;
1208     dlen -= pos;
1209     if (dlen <= 0) {
1210         strcpy(dst, src);
1211         return;
1212     }
1213     slen = strlen(src);
1214     do {
1215         dst[slen] = *dst;
1216         --dst;
1217     } while (--dlen >= 0);
1218     dst++;
1219     while (--slen >= 0)
1220         *dst++ = *src++;
1221 }
1222
1223
1224
1225
1226 /* File functions */
1227
1228 /* Peek at next character of input stream; return EOF at end-of-file. */
1229
1230 int P_peek(f)
1231 FILE *f;
1232 {
1233     int ch;
1234
1235     ch = getc(f);
1236     if (ch == EOF)
1237         return EOF;
1238     ungetc(ch, f);
1239     return (ch == '\n') ? ' ' : ch;
1240 }
1241
1242
1243 /* Check if at end of file, using Pascal "eof" semantics.  End-of-file for
1244    stdin is broken; remove the special case for it to be broken in a
1245    different way. */
1246
1247 int P_eof(f)
1248 FILE *f;
1249 {
1250     register int ch;
1251
1252     if (feof(f))
1253         return 1;
1254     if (f == stdin)
1255         return 0;    /* not safe to look-ahead on the keyboard! */
1256     ch = getc(f);
1257     if (ch == EOF)
1258         return 1;
1259     ungetc(ch, f);
1260     return 0;
1261 }
1262
1263
1264 /* Check if at end of line (or end of entire file). */
1265
1266 int P_eoln(f)
1267 FILE *f;
1268 {
1269     register int ch;
1270
1271     ch = getc(f);
1272     if (ch == EOF)
1273         return 1;
1274     ungetc(ch, f);
1275     return (ch == '\n');
1276 }
1277
1278
1279 /* Read a packed array of characters from a file. */
1280
1281 Void P_readpaoc(f, s, len)
1282 FILE *f;
1283 char *s;
1284 int len;
1285 {
1286     int ch;
1287
1288     for (;;) {
1289         if (len <= 0)
1290             return;
1291         ch = getc(f);
1292         if (ch == EOF || ch == '\n')
1293             break;
1294         *s++ = ch;
1295         --len;
1296     }
1297     while (--len >= 0)
1298         *s++ = ' ';
1299     if (ch != EOF)
1300         ungetc(ch, f);
1301 }
1302
1303 Void P_readlnpaoc(f, s, len)
1304 FILE *f;
1305 char *s;
1306 int len;
1307 {
1308     int ch;
1309
1310     for (;;) {
1311         ch = getc(f);
1312         if (ch == EOF || ch == '\n')
1313             break;
1314         if (len > 0) {
1315             *s++ = ch;
1316             --len;
1317         }
1318     }
1319     while (--len >= 0)
1320         *s++ = ' ';
1321 }
1322
1323
1324 /* Compute maximum legal "seek" index in file (0-based). */
1325
1326 long P_maxpos(f)
1327 FILE *f;
1328 {
1329     long savepos = ftell(f);
1330     long val;
1331
1332     if (fseek(f, 0L, SEEK_END))
1333         return -1;
1334     val = ftell(f);
1335     if (fseek(f, savepos, SEEK_SET))
1336         return -1;
1337     return val;
1338 }
1339
1340
1341 /* Use packed array of char for a file name. */
1342
1343 Char *P_trimname(fn, len)
1344 register Char *fn;
1345 register int len;
1346 {
1347     static Char fnbuf[256];
1348     register Char *cp = fnbuf;
1349
1350     while (--len >= 0 && *fn && !isspace(*fn))
1351         *cp++ = *fn++;
1352     *cp = 0;
1353     return fnbuf;
1354 }
1355
1356
1357
1358
1359
1360
1361 /* Sets are stored as an array of longs.  S[0] is the size of the set;
1362    S[N] is the N'th 32-bit chunk of the set.  S[0] equals the maximum
1363    I such that S[I] is nonzero.  S[0] is zero for an empty set.  Within
1364    each long, bits are packed from lsb to msb.  The first bit of the
1365    set is the element with ordinal value 0.  (Thus, for a "set of 5..99",
1366    the lowest five bits of the first long are unused and always zero.) */
1367
1368 /* (Sets with 32 or fewer elements are normally stored as plain longs.) */
1369
1370 long *P_setunion(d, s1, s2)         /* d := s1 + s2 */
1371 register long *d, *s1, *s2;
1372 {
1373     long *dbase = d++;
1374     register int sz1 = *s1++, sz2 = *s2++;
1375     while (sz1 > 0 && sz2 > 0) {
1376         *d++ = *s1++ | *s2++;
1377         sz1--, sz2--;
1378     }
1379     while (--sz1 >= 0)
1380         *d++ = *s1++;
1381     while (--sz2 >= 0)
1382         *d++ = *s2++;
1383     *dbase = d - dbase - 1;
1384     return dbase;
1385 }
1386
1387
1388 long *P_setint(d, s1, s2)           /* d := s1 * s2 */
1389 register long *d, *s1, *s2;
1390 {
1391     long *dbase = d++;
1392     register int sz1 = *s1++, sz2 = *s2++;
1393     while (--sz1 >= 0 && --sz2 >= 0)
1394         *d++ = *s1++ & *s2++;
1395     while (--d > dbase && !*d) ;
1396     *dbase = d - dbase;
1397     return dbase;
1398 }
1399
1400
1401 long *P_setdiff(d, s1, s2)          /* d := s1 - s2 */
1402 register long *d, *s1, *s2;
1403 {
1404     long *dbase = d++;
1405     register int sz1 = *s1++, sz2 = *s2++;
1406     while (--sz1 >= 0 && --sz2 >= 0)
1407         *d++ = *s1++ & ~*s2++;
1408     if (sz1 >= 0) {
1409         while (sz1-- >= 0)
1410             *d++ = *s1++;
1411     }
1412     while (--d > dbase && !*d) ;
1413     *dbase = d - dbase;
1414     return dbase;
1415 }
1416
1417
1418 long *P_setxor(d, s1, s2)         /* d := s1 / s2 */
1419 register long *d, *s1, *s2;
1420 {
1421     long *dbase = d++;
1422     register int sz1 = *s1++, sz2 = *s2++;
1423     while (sz1 > 0 && sz2 > 0) {
1424         *d++ = *s1++ ^ *s2++;
1425         sz1--, sz2--;
1426     }
1427     while (--sz1 >= 0)
1428         *d++ = *s1++;
1429     while (--sz2 >= 0)
1430         *d++ = *s2++;
1431     while (--d > dbase && !*d) ;
1432     *dbase = d - dbase;
1433     return dbase;
1434 }
1435
1436
1437 int P_inset(val, s)                 /* val IN s */
1438 register unsigned val;
1439 register long *s;
1440 {
1441     register int bit;
1442     bit = val % SETBITS;
1443     val /= SETBITS;
1444     if (val < *s++ && ((1<<bit) & s[val]))
1445         return 1;
1446     return 0;
1447 }
1448
1449
1450 long *P_addset(s, val)              /* s := s + [val] */
1451 register long *s;
1452 register unsigned val;
1453 {
1454     register long *sbase = s;
1455     register int bit, size;
1456     bit = val % SETBITS;
1457     val /= SETBITS;
1458     size = *s;
1459     if (++val > size) {
1460         s += size;
1461         while (val > size)
1462             *++s = 0, size++;
1463         *sbase = size;
1464     } else
1465         s += val;
1466     *s |= 1<<bit;
1467     return sbase;
1468 }
1469
1470
1471 long *P_addsetr(s, v1, v2)              /* s := s + [v1..v2] */
1472 register long *s;
1473 register unsigned v1, v2;
1474 {
1475     register long *sbase = s;
1476     register int b1, b2, size;
1477     if ((int)v1 > (int)v2)
1478         return sbase;
1479     b1 = v1 % SETBITS;
1480     v1 /= SETBITS;
1481     b2 = v2 % SETBITS;
1482     v2 /= SETBITS;
1483     size = *s;
1484     v1++;
1485     if (++v2 > size) {
1486         while (v2 > size)
1487             s[++size] = 0;
1488         s[v2] = 0;
1489         *s = v2;
1490     }
1491     s += v1;
1492     if (v1 == v2) {
1493         *s |= (~((-2)<<(b2-b1))) << b1;
1494     } else {
1495         *s++ |= (-1) << b1;
1496         while (++v1 < v2)
1497             *s++ = -1;
1498         *s |= ~((-2) << b2);
1499     }
1500     return sbase;
1501 }
1502
1503
1504 long *P_remset(s, val)              /* s := s - [val] */
1505 register long *s;
1506 register unsigned val;
1507 {
1508     register int bit;
1509     bit = val % SETBITS;
1510     val /= SETBITS;
1511     if (++val <= *s) {
1512         if (!(s[val] &= ~(1<<bit)))
1513             while (*s && !s[*s])
1514                 (*s)--;
1515     }
1516     return s;
1517 }
1518
1519
1520 int P_setequal(s1, s2)              /* s1 = s2 */
1521 register long *s1, *s2;
1522 {
1523     register int size = *s1++;
1524     if (*s2++ != size)
1525         return 0;
1526     while (--size >= 0) {
1527         if (*s1++ != *s2++)
1528             return 0;
1529     }
1530     return 1;
1531 }
1532
1533
1534 int P_subset(s1, s2)                /* s1 <= s2 */
1535 register long *s1, *s2;
1536 {
1537     register int sz1 = *s1++, sz2 = *s2++;
1538     if (sz1 > sz2)
1539         return 0;
1540     while (--sz1 >= 0) {
1541         if (*s1++ & ~*s2++)
1542             return 0;
1543     }
1544     return 1;
1545 }
1546
1547
1548 long *P_setcpy(d, s)                /* d := s */
1549 register long *d, *s;
1550 {
1551     register long *save_d = d;
1552
1553 #ifdef SETCPY_MEMCPY
1554     memcpy(d, s, (*s + 1) * sizeof(long));
1555 #else
1556     register int i = *s + 1;
1557     while (--i >= 0)
1558         *d++ = *s++;
1559 #endif
1560     return save_d;
1561 }
1562
1563
1564 /* s is a "smallset", i.e., a 32-bit or less set stored
1565    directly in a long. */
1566
1567 long *P_expset(d, s)                /* d := s */
1568 register long *d;
1569 register long s;
1570 {
1571     if (s) {
1572         d[1] = s;
1573         *d = 1;
1574     } else
1575         *d = 0;
1576     return d;
1577 }
1578
1579
1580 long P_packset(s)                   /* convert s to a small-set */
1581 register long *s;
1582 {
1583     if (*s++)
1584         return *s;
1585     else
1586         return 0;
1587 }
1588
1589
1590
1591
1592
1593 /* Oregon Software Pascal extensions, courtesy of William Bader */
1594
1595 int P_getcmdline(l, h, line)
1596 int l, h;
1597 Char *line;
1598 {
1599     int i, len;
1600     char *s;
1601
1602     h = h - l + 1;
1603     len = 0;
1604     for(i = 1; i < P_argc; i++) {
1605         s = P_argv[i];
1606         while (*s) {
1607             if (len >= h) return len;
1608             line[len++] = *s++;
1609         }
1610         if (len >= h) return len;
1611         line[len++] = ' ';
1612     }
1613     return len;
1614 }
1615
1616 Void TimeStamp(Day, Month, Year, Hour, Min, Sec)
1617 int *Day, *Month, *Year, *Hour, *Min, *Sec;
1618 {
1619 #ifndef NO_TIME
1620     struct tm *tm;
1621     long clock;
1622
1623     time(&clock);
1624     tm = localtime(&clock);
1625     *Day = tm->tm_mday;
1626     *Month = tm->tm_mon + 1;            /* Jan = 0 */
1627     *Year = tm->tm_year;
1628     if (*Year < 1900)
1629         *Year += 1900;     /* year since 1900 */
1630     *Hour = tm->tm_hour;
1631     *Min = tm->tm_min;
1632     *Sec = tm->tm_sec;
1633 #endif
1634 }
1635
1636
1637
1638
1639 /* SUN Berkeley Pascal extensions */
1640
1641 Void P_sun_argv(s, len, n)
1642 register char *s;
1643 register int len, n;
1644 {
1645     register char *cp;
1646
1647     if ((unsigned)n < P_argc)
1648       cp = P_argv[n];
1649     else
1650       cp = "";
1651     while (*cp && --len >= 0)
1652       *s++ = *cp++;
1653     while (--len >= 0)
1654       *s++ = ' ';
1655 }
1656
1657
1658
1659
1660 int _OutMem()
1661 {
1662     return _Escape(-2);
1663 }
1664
1665 int _CaseCheck()
1666 {
1667     return _Escape(-9);
1668 }
1669
1670 int _NilCheck()
1671 {
1672     return _Escape(-3);
1673 }
1674
1675
1676
1677
1678
1679 /* The following is suitable for the HP Pascal operating system.
1680    It might want to be revised when emulating another system. */
1681
1682 char *_ShowEscape(buf, code, ior, prefix)
1683 char *buf, *prefix;
1684 int code, ior;
1685 {
1686     char *bufp;
1687
1688     if (prefix && *prefix) {
1689         strcpy(buf, prefix);
1690         strcat(buf, ": ");
1691         bufp = buf + strlen(buf);
1692     } else {
1693         bufp = buf;
1694     }
1695     if (code == -10) {
1696         sprintf(bufp, "Pascal system I/O error %d", ior);
1697         switch (ior) {
1698             case 3:
1699                 strcat(buf, " (illegal I/O request)");
1700                 break;
1701             case 7:
1702                 strcat(buf, " (bad file name)");
1703                 break;
1704             case FileNotFound:   /*10*/
1705                 strcat(buf, " (file not found)");
1706                 break;
1707             case FileNotOpen:    /*13*/
1708                 strcat(buf, " (file not open)");
1709                 break;
1710             case BadInputFormat: /*14*/
1711                 strcat(buf, " (bad input format)");
1712                 break;
1713             case 24:
1714                 strcat(buf, " (not open for reading)");
1715                 break;
1716             case 25:
1717                 strcat(buf, " (not open for writing)");
1718                 break;
1719             case 26:
1720                 strcat(buf, " (not open for direct access)");
1721                 break;
1722             case 28:
1723                 strcat(buf, " (string subscript out of range)");
1724                 break;
1725             case EndOfFile:      /*30*/
1726                 strcat(buf, " (end-of-file)");
1727                 break;
1728             case FileWriteError: /*38*/
1729                 strcat(buf, " (file write error)");
1730                 break;
1731         }
1732     } else {
1733         sprintf(bufp, "Pascal system error %d", code);
1734         switch (code) {
1735             case -2:
1736                 strcat(buf, " (out of memory)");
1737                 break;
1738             case -3:
1739                 strcat(buf, " (reference to NIL pointer)");
1740                 break;
1741             case -4:
1742                 strcat(buf, " (integer overflow)");
1743                 break;
1744             case -5:
1745                 strcat(buf, " (divide by zero)");
1746                 break;
1747             case -6:
1748                 strcat(buf, " (real math overflow)");
1749                 break;
1750             case -8:
1751                 strcat(buf, " (value range error)");
1752                 break;
1753             case -9:
1754                 strcat(buf, " (CASE value range error)");
1755                 break;
1756             case -12:
1757                 strcat(buf, " (bus error)");
1758                 break;
1759             case -20:
1760                 strcat(buf, " (stopped by user)");
1761                 break;
1762         }
1763     }
1764     return buf;
1765 }
1766
1767
1768 int _Escape(code)
1769 int code;
1770 {
1771     char buf[100];
1772
1773     P_escapecode = code;
1774     if (__top_jb) {
1775         __p2c_jmp_buf *jb = __top_jb;
1776         __top_jb = jb->next;
1777         longjmp(jb->jbuf, 1);
1778     }
1779     if (code == 0)
1780         exit(EXIT_SUCCESS);
1781     if (code == -1)
1782         exit(EXIT_FAILURE);
1783     fprintf(stderr, "%s\n", _ShowEscape(buf, P_escapecode, P_ioresult, ""));
1784     exit(EXIT_FAILURE);
1785 }
1786
1787 int _EscIO(code)
1788 int code;
1789 {
1790     P_ioresult = code;
1791     return _Escape(-10);
1792 }
1793
1794
1795
1796
1797 /* End. */
1798
1799
1800
1801
1802
1803 /************************************************************
1804 *                                                           *
1805 *          date.c                                           *
1806 *                                                           *
1807 ************************************************************/
1808
1809
1810
1811 /*#include <time.h>
1812 #include <stdio.h>*/
1813 static  char *mon[]={"JAN","FEB","MAR","APR","MAY","JUN",
1814                 "JUL","AUG","SEP","OCT","NOV","DEC"};
1815 /* PROCEDURE DATE(VAR DATESTRING:PACKED ARRAY[1..11] OF CHAR);EXTERN; */
1816 /* activate DATE by removing comment brackets if necessary */
1817 /***/
1818
1819 void Date(string)
1820 char* string;
1821 {
1822   time_t tt;
1823   struct tm *t;
1824   time(&tt);
1825   t=localtime(&tt);
1826   sprintf(string,"%d-%s-19%d\n",t->tm_mday,mon[t->tm_mon],t->tm_year);
1827 }
1828
1829
1830 /* p2c: dssp.p, line 295:
1831  * Note: Unexpected name "tapein" in program header [262] */
1832 /* p2c: dssp.p, line 295:
1833  * Note: Unexpected name "tapeout" in program header [262] */
1834
1835
1836 /*--------------------------------------------------------------------*/
1837 /* PROGRAM FATAL ERROR EXIT LABEL */
1838 /*******************  MATHEMATICAL CONSTANTS  **************************
1839  YVERTEX, - ARE Y,Z-COMPONENTS OF THE FIRST ICOSAHEDRON VERTEX. THE
1840  ZVERTEX    X-COMPONENT IS 0.
1841  EPS      - NUMERICAL TOLERANCE
1842   --------------------------------------------------------------------*/
1843
1844 #define PIHALF          1.570796
1845 #define PI              3.141593
1846 #define TWOPI           6.283185
1847 #define FOURPI          12.56637
1848 #define RADIAN          57.29578
1849 #define YVERTEX         0.8506508
1850 #define ZVERTEX         0.5257311
1851 #define EPS             0.00001
1852 /***/
1853 /***************  ARRAY DIMENSIONING CONSTANTS  ***********************
1854  NMAX     - MAXIMUM NUMBER OF AMINOACID RESIDUES IN ARRAY CHAIN
1855  MAXATOM  - MAXIMUM NUMBER OF SIDECHAIN ATOMS IN ARRAY SIDECHAIN
1856  MAXBRIDGE- MAXIMUM NUMBER OF BRIDGES IN ARRAY BRIDGETABLE
1857  NFACE,   - NUMBER OF FACES OF POLYHEDRON. THE COORDINATES OF THE CENTRE
1858  ORDER      OF EACH TRIANGULAR FACE ARE STORED IN ARRAY P, THE AREA
1859              IS STORED IN ARRAY WP IN PROCEDURE FLAGACCESS. NFACE MUST
1860              BE OF THE FORM NFACE=20*(4**ORDER), ORDER=0,1,2,...
1861              THE ACCURACY OF THE SOLVENT ACCESSIBLE SURFACE OF EACH
1862              AMINOACID RESIDUE IS ONE ANGSTROM**2 FOR ORDER=2,NFACE=320.
1863  MAXPACK  - MAXIMUM NUMBER OF PROTEIN ATOMS WHICH CAN INTRUDE INTO
1864              SOLVENT AROUND ANY GIVEN TEST ATOM. THE COORDINATES OF
1865              THESE ATOMS ARE STORED IN ARRAY X, THEIR RADII IN ARRAY RX
1866              IN PROCEDURE SURFACE.
1867  MAXHIST  - NUMBER OF SLOTS IN ARRAYS HELIXHIST AND BETAHIST USED FOR
1868              LENGTH STATISTICS OF SECONDARY STRUCTURE.
1869  MAXSS    - MAXIMUM NUMBER OF SSBOND RECORDS ON INPUT FILE. THE
1870              DISULFIDE BOND ARE SAVED IN ARRAY SSBONDS.
1871   --------------------------------------------------------------------*/
1872
1873 #define NMAX            6000
1874
1875 #define MAXATOM         40000L
1876
1877 #define MAXBRIDGE       300
1878 #define NFACE           320
1879 #define ORDER           2
1880 #define MAXPACK         200
1881 #define MAXHIST         30
1882 #define MAXSS           100
1883 /***/
1884 /*********************  PHYSICAL CONSTANTS   **************************
1885  RN       - RADIUS OF PEPTIDE NITROGEN ATOM
1886  RCA      - RADIUS OF PEPTIDE ALPHA-CARBON ATOM
1887  RC       - RADIUS OF PEPTIDE C'-CARBON ATOM
1888  RO       - RADIUS OF PEPTIDE OXYGEN ATOM
1889  RSIDEATOM- RADIUS OF SIDECHAIN ATOM
1890  RWATER   - RADIUS OF WATER MOLECULE
1891  SSDIST   - MAXIMUM ALLOWED DISTANCE OF DISULFIDE BRIDGE
1892  BREAKDIST- MAXIMUM ALLOWED PEPTIDE BOND LENGTH. IF DISTANCE IS
1893              GREATER A POLYPEPTIDE CHAIN INTERRUPTION IS ASSUMED.
1894  RESRAD   - MAXIMUM RADIUS OF A SPHERE AROUND C-ALPHA CONTAINING
1895              ALL ATOMS OF A RESIDUE
1896  CADIST   - MINIMUM DISTANCE BETWEEN ALPHA-CARBON ATOMS SUCH THAT NO
1897              BACKBONE HYDROGEN BONDS CAN BE FORMED
1898  DIST     - SMALLEST ALLOWED DISTANCE BETWEEN ANY ATOMS
1899  MAXDIST  - LARGEST ALLOWED DISTANCE BETWEEN SIDECHAIN ATOM AND C-ALPHA
1900              WITHIN A RESIDUE
1901  Q        - COUPLING CONSTANT FOR ELECTROSTATIC ENERGY
1902                     Q=-332*0.42*0.2*1000.0
1903  HBLOW    - LOWEST ALLOWED  ENERGY OF A HYDROGEN BOND IN CAL/MOL
1904  HBHIGH   - HIGHEST ALLOWED ENERGY OF A HYDROGEN BOND IN CAL/MOL
1905   --------------------------------------------------------------------*/
1906
1907 #define RN              1.65
1908 #define RCA             1.87
1909 #define RC              1.76
1910 #define RO              1.4
1911 #define RSIDEATOM       1.8
1912 #define RWATER          1.4
1913 #define SSDIST          3.0
1914 #define BREAKDIST       2.5
1915 #define RESRAD          10.0
1916 #define CADIST          9.0
1917 #define DIST            0.5
1918 #define MAXDIST         10.0
1919 #define Q               (-27888.0)
1920
1921 #define HBLOW           (-9900)
1922 #define HBHIGH          (-500)
1923
1924
1925 /***/
1926 /***************** GLOBAL DATA TYPE DEFINITIONS ***********************/
1927
1928 typedef double vector[3];
1929 typedef Char char4[4];
1930 typedef Char char6[6];
1931 typedef enum {
1932   parallel, antiparallel, nobridge
1933 } bridgetyp;
1934 typedef enum {
1935   symbol, turn3, turn4, turn5, bend, chirality, beta1, beta2
1936 } structure;
1937
1938 typedef struct hydrogenbond {
1939   long residue, energy;
1940 } hydrogenbond;
1941
1942 typedef hydrogenbond bonds[2];
1943
1944 typedef struct backbone {
1945   char6 aaident;
1946   Char sheetlabel, aa;
1947   char4 threelettercode;
1948   Char ss[(long)beta2 - (long)symbol + 1];
1949   long partner[(long)beta2 - (long)beta1 + 1];
1950   long access;
1951   double alpha, kappa;
1952   bonds acceptor, donor;
1953   vector h, n, ca, c, o;
1954   long atompointer, nsideatoms;
1955 } backbone;
1956
1957 typedef struct bridge {
1958   Char sheetname, laddername;
1959   bridgetyp btyp;
1960   long linkset[MAXBRIDGE / 32 + 2];
1961   long ib, ie, jb, je, from, towards;
1962 } bridge;
1963
1964
1965 static int bVerbose;
1966 Static int silentFlag;
1967 Static long nss, nssintra, nssinter, lchain, nbridge;
1968 Static char6 ssbonds[MAXSS][2];
1969 Static backbone chain[NMAX + 1];
1970 FILE *tapein, *tapeout;
1971 Static vector sidechain[MAXATOM];
1972 Static bridge bridgetable[MAXBRIDGE];
1973
1974 Static Void VecCopy(dest,source)
1975 double* dest;
1976 double* source;
1977 {
1978   dest[0]=source[0];
1979   dest[1]=source[1];
1980   dest[2]=source[2];
1981 }
1982
1983 Static Void StrCopy(dest,source,n)
1984 char* dest;
1985 char* source;
1986 int n;
1987 {
1988   int i;
1989   for(i=0;i<n;i++)
1990     dest[i]=source[i];
1991 }
1992
1993 Static double Atan2(y, x)
1994 double y, x;
1995 {
1996   double z;
1997
1998   if (x != 0.0)
1999     z = atan(y / x);
2000   else if (y > 0.0)
2001     z = PIHALF;
2002   else if (y < 0.0)
2003     z = -PIHALF;
2004   else
2005     z = TWOPI;
2006   if (x >= 0.0)
2007     return z;
2008   if (y > 0.0)
2009     z += PI;
2010   else
2011     z -= PI;
2012   return z;
2013 }  /* Atan2 */
2014
2015
2016 /***/
2017
2018 Static Void Diff(x, y, z)
2019 double *x, *y, *z;
2020 {
2021   z[0] = x[0] - y[0];
2022   z[1] = x[1] - y[1];
2023   z[2] = x[2] - y[2];
2024 }  /* Diff */
2025
2026
2027 /***/
2028
2029 Static double Dot(x, y)
2030 double *x, *y;
2031 {
2032   return (x[0] * y[0] + x[1] * y[1] + x[2] * y[2]);
2033 }  /* Dot */
2034
2035
2036 /***/
2037
2038 Static Void Cross(x, y, z)
2039 double *x, *y, *z;
2040 {
2041   z[0] = x[1] * y[2] - y[1] * x[2];
2042   z[1] = x[2] * y[0] - y[2] * x[0];
2043   z[2] = x[0] * y[1] - y[0] * x[1];
2044 }  /* Cross */
2045
2046
2047 /***/
2048
2049 Static Void Norm(x, xnorm)
2050 double *x;
2051 double *xnorm;
2052 {
2053   /* RETURNS INPUT VECTOR X NORMALIZED TO UNIT LENGTH.
2054      XNORM IS THE ORIGINAL LENGTH OF X.                         */
2055   double TEMP, TEMP1, TEMP2;
2056
2057   TEMP = x[0];
2058   TEMP1 = x[1];
2059   TEMP2 = x[2];
2060   *xnorm = TEMP * TEMP + TEMP1 * TEMP1 + TEMP2 * TEMP2;
2061   if (*xnorm <= 0.0)
2062     return;
2063   *xnorm = sqrt(*xnorm);
2064   x[0] /= *xnorm;
2065   x[1] /= *xnorm;
2066   x[2] /= *xnorm;
2067 }  /* Norm */
2068
2069
2070 /***/
2071
2072 Static double Dihedralangle(v1, v2, v3, v4)
2073 double *v1, *v2, *v3, *v4;
2074 {
2075   /*CALCULATES TORSION ANGLE OF A SET OF 4 ATOMS V1-V2-V3-V4.
2076     DIHEDRALANGLE IS THE ANGLE BETWEEN THE PROJECTION OF
2077     V1-V2 AND THE PROJECTION OF V4-V3 ONTO A PLANE NORMAL TO
2078     BOND V2-V3.*/
2079   /***/
2080   double Result, u, v;
2081   vector v12, v43, x, y, z, p;
2082
2083   Diff(v1, v2, v12);
2084   Diff(v4, v3, v43);
2085   Diff(v2, v3, z);
2086   Cross(z, v12, p);
2087   Cross(z, v43, x);
2088   Cross(z, x, y);
2089   u = Dot(x, x);
2090   v = Dot(y, y);
2091   Result = 360.0;
2092   if (u <= 0.0 || v <= 0.0)
2093     return Result;
2094   u = Dot(p, x) / sqrt(u);
2095   v = Dot(p, y) / sqrt(v);
2096   if (u != 0.0 || v != 0.0)
2097     return (Atan2(v, u) * RADIAN);
2098   return Result;
2099 }  /* Dihedralangle */
2100
2101
2102 /***/
2103
2104 Static double Cosangle(v1, v2, v3, v4)
2105 double *v1, *v2, *v3, *v4;
2106 {
2107   vector u, v;
2108   double x;
2109
2110   Diff(v1, v2, u);
2111   Diff(v3, v4, v);
2112   x = Dot(u, u) * Dot(v, v);
2113   if (x > 0.0)
2114     return (Dot(u, v) / sqrt(x));
2115   else
2116     return 0.0;
2117 }  /* Cosangle */
2118
2119
2120 /***/
2121
2122 Static double Distance(u, v)
2123 double *u, *v;
2124 {
2125   double TEMP, TEMP1, TEMP2;
2126
2127   TEMP = u[0] - v[0];
2128   TEMP1 = u[1] - v[1];
2129   TEMP2 = u[2] - v[2];
2130   return sqrt(TEMP * TEMP + TEMP1 * TEMP1 + TEMP2 * TEMP2);
2131 }  /* Distance */
2132
2133
2134 /***/
2135
2136 Static double Distsq(u, v)
2137 double *u, *v;
2138 {
2139   double TEMP, TEMP1, TEMP2;
2140
2141   TEMP = u[0] - v[0];
2142   TEMP1 = u[1] - v[1];
2143   TEMP2 = u[2] - v[2];
2144   return (TEMP * TEMP + TEMP1 * TEMP1 + TEMP2 * TEMP2);
2145 }  /* Distsq */
2146
2147
2148 /*--------------------------------------------------------------------*/
2149
2150 Static boolean Nochainbreak(i, j)
2151 long i, j;
2152 {
2153   long k;
2154   boolean test;
2155
2156   test = (i >= 1 && j <= NMAX && i <= j);
2157   k = i;
2158   while (test && k <= j) {
2159     if (chain[k].aa == '!')
2160       test = false;
2161     else
2162       k++;
2163   }
2164   return test;
2165 }  /* Nochainbreak */
2166
2167
2168 /***/
2169 /*--------------------------------------------------------------------*/
2170
2171 Static Void Writeresidue(res)
2172 backbone res;
2173 {
2174   long i;
2175
2176   for (i = 0; i <= 3; i++)
2177     putchar(res.threelettercode[i]);
2178   for (i = 0; i <= 5; i++)
2179     putchar(res.aaident[i]);
2180 }  /* Writeresidue */
2181
2182
2183 #define MAXSIDEATOMS    20
2184
2185
2186 typedef enum {
2187   headercard, compndcard, sourcecard, authorcard, ssbondcard, atomcard,
2188   tercard, endcard, othercard
2189 } cardtype;
2190 /***/
2191
2192 typedef struct cardcontents {
2193   cardtype art;
2194   union {
2195     Char z[128];
2196     char6 r[2];
2197     struct {
2198       char4 atomname, aaname;
2199       Char altloc, residuename;
2200       char6 reseqnum;
2201       vector coordinates;
2202     } U5;
2203     Char ch;
2204   } UU;
2205 } cardcontents;   /* CARDCONTENTS TYPE DEFINITION */
2206
2207 /***/
2208
2209
2210 Static jmp_buf _JL99;
2211
2212 /* Local variables for Inputcoordinates: */
2213 struct LOC_Inputcoordinates {
2214   long *lchain, latom, hatoms;
2215   boolean nmissing, camissing, cmissing, omissing, corelimit;
2216   vector sidecoordinates[MAXSIDEATOMS];
2217   double dco;
2218   char4 sideatomnames[MAXSIDEATOMS];
2219   backbone reszero, resinfo;
2220 } ;
2221
2222 /***/
2223
2224 Local Char Onelettercode(aaa, LINK)
2225 Char *aaa;
2226 struct LOC_Inputcoordinates *LINK;
2227 {
2228   Char aasymbol[50];
2229   Char aminoacid[150];
2230   Char string[5][30];
2231   long i, l, k;
2232   Char a;
2233
2234   StrCopy(aasymbol, "ARNDCEQGHILKMFPSTWYVBZXXXXXXXXXXXXXXXX--CCCCIPPPW-", 50L);
2235   StrCopy(string[0], "ALAARGASNASPCYSGLUGLNGLYHISILE", 30L);
2236   StrCopy(string[1], "LEULYSMETPHEPROSERTHRTRPTYRVAL", 30L);
2237   StrCopy(string[2], "ASXGLXACDALBALIABUAROBASBETHSE", 30L);
2238   StrCopy(string[3], "HYPHYLORNPCASARTAUTHYUNKACEFOR", 30L);
2239   StrCopy(string[4], "CYHCSHCSSCYXILUPRZPR0CPRTRYHOH", 30L);
2240   l = 0;
2241   for (k = 0; k <= 4; k++) {
2242     for (i = 0; i <= 29; i++) {
2243       l++;
2244       aminoacid[l - 1] = string[k][i];
2245     }
2246   }
2247   a = '-';
2248   i = 1;
2249   k = 1;
2250   while (k < 51 && a == '-') {
2251     if (aminoacid[i - 1] == aaa[0]) {
2252       if (aminoacid[i] == aaa[1]) {
2253         if (aminoacid[i + 1] == aaa[2])
2254           a = aasymbol[k - 1];
2255       }
2256     }
2257     i += 3;
2258     k++;
2259   }
2260   return a;
2261 }  /* Onelettercode */
2262
2263 /* Local variables for Checksideatoms: */
2264 struct LOC_Checksideatoms {
2265   struct LOC_Inputcoordinates *LINK;
2266 } ;
2267
2268 /***/
2269
2270 Local Void Checkdist(resinfo, LINK)
2271 backbone *resinfo;
2272 struct LOC_Checksideatoms *LINK;
2273 {
2274   long i, j, FORLIM;
2275
2276   i = 1;
2277   while (i <= resinfo->nsideatoms) {
2278     if (Distance(resinfo->ca, LINK->LINK->sidecoordinates[i - 1]) <= MAXDIST) {
2279       i++;
2280       continue;
2281     }
2282     if (bVerbose) {
2283       printf(" !!! RESIDUE ");
2284       Writeresidue(*resinfo);
2285       printf(" HAS ILLEGAL SIDECHAIN ATOM NAMED ");
2286       for (j = 0; j <= 3; j++)
2287         putchar(LINK->LINK->sideatomnames[i - 1][j]);
2288       printf(".\n");
2289       printf("     THIS ATOM WILL BE IGNORED !!!\n\n");
2290     }
2291     FORLIM = resinfo->nsideatoms;
2292     for (j = i + 1; j <= FORLIM; j++) {
2293       StrCopy(LINK->LINK->sideatomnames[j - 2],
2294               LINK->LINK->sideatomnames[j - 1], sizeof(char4));
2295       VecCopy(LINK->LINK->sidecoordinates[j - 2],
2296               LINK->LINK->sidecoordinates[j - 1]);
2297     }
2298     resinfo->nsideatoms--;
2299   }
2300 }  /* Checkdist */
2301
2302 /***/
2303
2304 Local Void Checksideatoms(resinfo, LINK)
2305 backbone *resinfo;
2306 struct LOC_Inputcoordinates *LINK;
2307 {
2308   struct LOC_Checksideatoms V;
2309   long i, j;
2310   Char c;
2311
2312   /***/
2313
2314   V.LINK = LINK;
2315   Checkdist(resinfo, &V);
2316   i = -1;
2317   c = resinfo->aa;
2318   if (c == 'G')
2319     i = 0;
2320   if (c == 'A')
2321     i = 1;
2322   if (c == 'S' || c == 'C')
2323     i = 2;
2324   if (c == 'P' || c == 'T' || c == 'V')
2325     i = 3;
2326   if (c == 'B' || c == 'M' || c == 'L' || c == 'I' || c == 'D' || c == 'N')
2327     i = 4;
2328   if (c == 'Z' || c == 'K' || c == 'Q' || c == 'E')
2329     i = 5;
2330   if (c == 'H')
2331     i = 6;
2332   if (c == 'F' || c == 'R')
2333     i = 7;
2334   if (c == 'Y')
2335     i = 8;
2336   if (c == 'W')
2337     i = 10;
2338   if ((resinfo->nsideatoms < i) && (bVerbose)) {
2339       printf(" !!! RESIDUE ");
2340       Writeresidue(*resinfo);
2341       printf(" HAS%3ld INSTEAD OF EXPECTED ", resinfo->nsideatoms);
2342       printf("%3ld SIDECHAIN ATOMS.\n", i);
2343       printf("     CALCULATED SOLVENT ACCESSIBILITY REFERS TO INCOMPLETE "
2344              "SIDECHAIN !!!\n\n");
2345     }
2346   if (i == -1 || resinfo->nsideatoms <= i)
2347     return;
2348   if (bVerbose) {
2349     printf(" !!! RESIDUE ");
2350     Writeresidue(*resinfo);
2351     printf(" HAS%3ld INSTEAD OF EXPECTED ", resinfo->nsideatoms);
2352     printf("%3ld SIDECHAIN ATOMS.\n", i);
2353     printf("     LAST SIDECHAIN ATOM NAME IS ");
2354     for (j = 0; j <= 3; j++)
2355       putchar(LINK->sideatomnames[resinfo->nsideatoms - 1][j]);
2356     printf("\n     CALCULATED SOLVENT ACCESSIBILITY INCLUDES EXTRA ATOMS !!!\n\n");
2357   }
2358 }  /* Checksideatoms */
2359
2360 /***/
2361
2362 Local Void Putresidue(LINK)
2363 struct LOC_Inputcoordinates *LINK;
2364 {
2365   /* insert residue into protein chain */
2366   long i;
2367   boolean complete;
2368   long FORLIM;
2369
2370   complete = !(LINK->nmissing || LINK->camissing || LINK->cmissing ||
2371                LINK->omissing);
2372   if (!complete &&
2373       strncmp(LINK->reszero.aaident, LINK->resinfo.aaident, sizeof(char6))
2374       && bVerbose) {
2375     printf(" !!! BACKBONE INCOMPLETE FOR RESIDUE ");
2376     Writeresidue(LINK->resinfo);
2377     printf("\n     RESIDUE WILL BE IGNORED !!!\n\n");
2378   }
2379   LINK->corelimit = (LINK->latom + LINK->resinfo.nsideatoms > MAXATOM ||
2380                      *LINK->lchain > NMAX - 2);
2381   if (complete && !LINK->corelimit) {
2382     Checksideatoms(&LINK->resinfo, LINK);
2383     VecCopy(LINK->resinfo.h, LINK->resinfo.n);
2384     if (Nochainbreak(*LINK->lchain, *LINK->lchain)) {
2385       if (Distance(chain[*LINK->lchain].c, LINK->resinfo.n) > BREAKDIST)
2386         /* keep ! at LCHAIN */
2387         /* CS Oct 1987 */
2388         if (bVerbose) {
2389           printf(" !!! EXCESSIVE C TO N DISTANCE ");
2390           printf("% .5E>% .5E\n",
2391                  Distance(chain[*LINK->lchain].c, LINK->resinfo.n), BREAKDIST);
2392           printf("     BEFORE RESIDUE ");
2393           Writeresidue(LINK->resinfo);
2394           printf(". CHAIN BREAK RESIDUE INSERTED !!!\n\n");
2395           (*LINK->lchain)++;
2396         }
2397     }
2398     if (Nochainbreak(*LINK->lchain, *LINK->lchain) && LINK->resinfo.aa != 'P') {
2399       LINK->dco = Distance(chain[*LINK->lchain].c, chain[*LINK->lchain].o);
2400       for (i = 0; i <= 2; i++)
2401         LINK->resinfo.h[i] = LINK->resinfo.n[i] +
2402             (chain[*LINK->lchain].c[i] - chain[*LINK->lchain].o[i]) / LINK->dco;
2403     }
2404     (*LINK->lchain)++;
2405     chain[*LINK->lchain] = LINK->resinfo;
2406     FORLIM = LINK->resinfo.nsideatoms;
2407     for (i = 0; i < FORLIM; i++)
2408       VecCopy(sidechain[LINK->latom + i], LINK->sidecoordinates[i]);
2409     LINK->latom += LINK->resinfo.nsideatoms;
2410   }
2411   if (Nochainbreak(*LINK->lchain, *LINK->lchain) && !complete)
2412     (*LINK->lchain)++;
2413   LINK->resinfo = LINK->reszero;
2414   LINK->nmissing = true;
2415   LINK->camissing = true;
2416   LINK->cmissing = true;
2417   LINK->omissing = true;
2418 }  /* Putresidue */
2419
2420 /***/
2421
2422 Local Void Getresidue(atomname, coordinates, LINK)
2423 Char *atomname;
2424 double *coordinates;
2425 struct LOC_Inputcoordinates *LINK;
2426 {
2427   boolean hydrogenatom;
2428
2429   hydrogenatom = ((atomname[0] == '9' || atomname[0] == '8' ||
2430                    atomname[0] == '7' || atomname[0] == '6' ||
2431                    atomname[0] == '5' || atomname[0] == '4' ||
2432                    atomname[0] == '3' || atomname[0] == '2' ||
2433                    atomname[0] == '1' || atomname[0] == '0' ||
2434                    atomname[0] == ' ') &&
2435                   (atomname[1] == 'D' || atomname[1] == 'H'));
2436   if (hydrogenatom) {
2437     LINK->hatoms++;
2438     return;
2439   }
2440   if (!strncmp(atomname, " N  ", sizeof(char4))) {
2441     LINK->nmissing = false;
2442     VecCopy(LINK->resinfo.n, coordinates);
2443     return;
2444   }
2445   if (!strncmp(atomname, " CA ", sizeof(char4))) {
2446     LINK->camissing = false;
2447     VecCopy(LINK->resinfo.ca, coordinates);
2448     return;
2449   }
2450   if (!strncmp(atomname, " C  ", sizeof(char4))) {
2451     LINK->cmissing = false;
2452     VecCopy(LINK->resinfo.c, coordinates);
2453     return;
2454   }
2455   if (!strncmp(atomname, " O  ", sizeof(char4))) {
2456     LINK->omissing = false;
2457     VecCopy(LINK->resinfo.o, coordinates);
2458     return;
2459   }
2460   if (LINK->resinfo.nsideatoms >= MAXSIDEATOMS)
2461     return;
2462   LINK->resinfo.nsideatoms++;
2463   VecCopy(LINK->sidecoordinates[LINK->resinfo.nsideatoms - 1], coordinates);
2464   StrCopy(LINK->sideatomnames[LINK->resinfo.nsideatoms - 1], atomname,
2465          sizeof(char4));
2466 }  /* Getresidue */
2467
2468 /***/
2469
2470 Local Void Readcard(cardinfo, LINK)
2471 cardcontents *cardinfo;
2472 struct LOC_Inputcoordinates *LINK;
2473 {
2474   Char c;
2475   long k, l, m;
2476   char6 key;
2477
2478   cardinfo->art = othercard;
2479   do {
2480     if (!P_eof(tapein)) {
2481       *key = getc(tapein);
2482       if (key[0] == '\n')
2483         key[0] = ' ';
2484     }
2485   } while (!(isupper(key[0]) | P_eof(tapein)));
2486   if (P_eof(tapein)) {
2487     cardinfo->art = endcard;
2488     return;
2489   }
2490   for (l = 1; l <= 5; l++) {
2491     if (!P_eoln(tapein)) {
2492       key[l] = getc(tapein);
2493       if (key[l] == '\n')
2494         key[l] = ' ';
2495     }
2496   }
2497   if (!strncmp(key, "HEADER", sizeof(char6)))
2498     cardinfo->art = headercard;
2499   if (!strncmp(key, "COMPND", sizeof(char6)))
2500     cardinfo->art = compndcard;
2501   if (!strncmp(key, "SOURCE", sizeof(char6)))
2502     cardinfo->art = sourcecard;
2503   if (!strncmp(key, "AUTHOR", sizeof(char6)))
2504     cardinfo->art = authorcard;
2505   if (!strncmp(key, "SSBOND", sizeof(char6)))
2506     cardinfo->art = ssbondcard;
2507   if (!strncmp(key, "ATOM  ", sizeof(char6)))
2508     cardinfo->art = atomcard;
2509   if (!strncmp(key, "TER   ", sizeof(char6)))
2510     cardinfo->art = tercard;
2511   if (!strncmp(key, "END   ", sizeof(char6)))
2512     cardinfo->art = endcard;
2513   switch (cardinfo->art) {
2514
2515   case headercard:
2516   case compndcard:
2517   case sourcecard:
2518   case authorcard:
2519     for (l = 0; l <= 5; l++)
2520       cardinfo->UU.z[l] = key[l];
2521     for (l = 6; l <= 126; l++)
2522       cardinfo->UU.z[l] = ' ';
2523     cardinfo->UU.z[127] = '.';
2524     if (cardinfo->art == headercard)
2525       m = 66;
2526     else
2527       m = 70;
2528     for (l = 6; l < m; l++) {
2529       if (!P_eoln(tapein)) {
2530         cardinfo->UU.z[l] = getc(tapein);
2531         if (cardinfo->UU.z[l] == '\n')
2532           cardinfo->UU.z[l] = ' ';
2533       }
2534     }
2535     break;
2536
2537   case ssbondcard:
2538     for (l = 7; l <= 8; l++) {
2539       c = getc(tapein);
2540       if (c == '\n')
2541         c = ' ';
2542     }
2543     for (k = 0; k <= 1; k++) {
2544       for (l = 1; l <= 7; l++) {
2545         c = getc(tapein);
2546         if (c == '\n')
2547           c = ' ';
2548       }
2549       cardinfo->UU.r[k][5] = getc(tapein);
2550       c = getc(tapein);
2551       if (cardinfo->UU.r[k][5] == '\n')
2552         cardinfo->UU.r[k][5] = ' ';
2553       if (c == '\n')
2554         c = ' ';
2555       /* minor modification suggested by Steven Sheriff */
2556       for (l = 0; l <= 3; l++) {
2557         cardinfo->UU.r[k][l] = getc(tapein);
2558         if (cardinfo->UU.r[k][l] == '\n')
2559           cardinfo->UU.r[k][l] = ' ';
2560       }
2561       if (P_eoln(tapein))
2562         cardinfo->UU.r[k][4] = ' ';
2563       else {
2564         cardinfo->UU.r[k][4] = getc(tapein);
2565         if (cardinfo->UU.r[k][4] == '\n')
2566           cardinfo->UU.r[k][4] = ' ';
2567       }
2568     }
2569     /* end minor modification suggested by Steven Sheriff */
2570     break;
2571
2572   case atomcard:
2573     for (l = 7; l <= 12; l++) {
2574       c = getc(tapein);
2575       if (c == '\n')
2576         c = ' ';
2577     }
2578     for (l = 0; l <= 3; l++) {
2579       cardinfo->UU.U5.atomname[l] = getc(tapein);
2580       if (cardinfo->UU.U5.atomname[l] == '\n')
2581         cardinfo->UU.U5.atomname[l] = ' ';
2582     }
2583     cardinfo->UU.U5.altloc = getc(tapein);
2584     if (cardinfo->UU.U5.altloc == '\n')
2585       cardinfo->UU.U5.altloc = ' ';
2586     for (l = 0; l <= 2; l++) {
2587       cardinfo->UU.U5.aaname[l] = getc(tapein);
2588       if (cardinfo->UU.U5.aaname[l] == '\n')
2589         cardinfo->UU.U5.aaname[l] = ' ';
2590     }
2591     cardinfo->UU.U5.aaname[3] = ' ';
2592     cardinfo->UU.U5.residuename = Onelettercode(cardinfo->UU.U5.aaname, LINK);
2593     c = getc(tapein);
2594     cardinfo->UU.U5.reseqnum[5] = getc(tapein);
2595     if (c == '\n')
2596       c = ' ';
2597     if (cardinfo->UU.U5.reseqnum[5] == '\n')
2598       cardinfo->UU.U5.reseqnum[5] = ' ';
2599     for (l = 0; l <= 4; l++) {
2600       cardinfo->UU.U5.reseqnum[l] = getc(tapein);
2601       if (cardinfo->UU.U5.reseqnum[l] == '\n')
2602         cardinfo->UU.U5.reseqnum[l] = ' ';
2603     }
2604     for (l = 0; l <= 2; l++)
2605       fscanf(tapein, "%lf", &cardinfo->UU.U5.coordinates[l]);
2606     break;
2607
2608   case tercard:
2609   case endcard:
2610   case othercard:
2611     /* blank case */
2612     break;
2613   }
2614   fscanf(tapein, "%*[^\n]");
2615   getc(tapein);
2616 }  /* Readcard */
2617
2618
2619 /***/
2620 /*--------------------------------------------------------------------*/
2621 /* SEE BROOKHAVEN PROTEIN DATA BANK ATOMIC COORDINATE ENTRY FORMAT
2622                                     OF DEC. 1981.
2623    -------------------------------------------------------------------*/
2624
2625 Static Void Inputcoordinates(lchain_)
2626 long *lchain_;
2627 {
2628   struct LOC_Inputcoordinates V;
2629   Char datestring[30];
2630   long i, j;
2631   boolean finish;
2632   structure s;
2633   cardtype ctype;
2634   cardcontents cardinfo;
2635   long cardhist[(long)othercard - (long)headercard + 1];
2636
2637   /***/
2638
2639   V.lchain = lchain_;
2640   nss = 0;
2641   V.latom = 0;
2642   V.hatoms = 0;
2643   for (j = 0; j <= 5; j++)
2644     V.reszero.aaident[j] = ' ';
2645   V.reszero.aa = '!';
2646   V.reszero.access = 0;
2647   StrCopy(V.reszero.threelettercode, "    ", sizeof(char4));
2648   for (s = symbol; (long)s <= (long)beta2; s = (structure)((long)s + 1))
2649     V.reszero.ss[(long)s - (long)symbol] = ' ';
2650   V.reszero.sheetlabel = ' ';
2651   V.reszero.partner[0] = 0;
2652   V.reszero.partner[(long)beta2 - (long)beta1] = 0;
2653   V.reszero.alpha = 360.0;
2654   V.reszero.kappa = 360.0;
2655   for (j = 0; j <= 1; j++) {
2656     V.reszero.acceptor[j].residue = 0;
2657     V.reszero.acceptor[j].energy = 0;
2658     V.reszero.donor[j].residue = 0;
2659     V.reszero.donor[j].energy = 0;
2660   }
2661   V.reszero.atompointer = 0;
2662   V.reszero.nsideatoms = 0;
2663   for (j = 0; j <= 2; j++) {
2664     V.reszero.h[j] = 0.0;
2665     V.reszero.n[j] = 0.0;
2666     V.reszero.ca[j] = 0.0;
2667     V.reszero.c[j] = 0.0;
2668     V.reszero.o[j] = 0.0;
2669   }
2670   for (i = 0; i <= NMAX; i++)
2671     chain[i] = V.reszero;
2672   Date(datestring);   /* DATE(DAY-MONTH-YEAR); */
2673   /* comment out this line if necessary */
2674   fprintf(tapeout, "**** SECONDARY STRUCTURE DEFINITION ");
2675   fprintf(tapeout, "BY THE PROGRAM DSSP, VERSION OCT. 1988 ****");
2676   fprintf(tapeout, " DATE=%.11s", datestring);
2677   for (i = 106; i <= 127; i++)
2678     putc(' ', tapeout);
2679   fprintf(tapeout, ".\n");
2680   fprintf(tapeout, "REFERENCE W. KABSCH AND C.SANDER, BIOPOLYMERS ");
2681   fprintf(tapeout, "22 (1983) 2577-2637");
2682   for (i = 66; i <= 127; i++)
2683     putc(' ', tapeout);
2684   fprintf(tapeout, ".\n");
2685   for (ctype = headercard;
2686        (long)ctype <= (long)othercard;
2687        ctype = (cardtype)((long)ctype + 1))
2688     cardhist[(long)ctype - (long)headercard] = 0;
2689   V.corelimit = false;
2690   finish = false;
2691   V.resinfo = V.reszero;
2692   V.nmissing = true;
2693   V.camissing = true;
2694   V.cmissing = true;
2695   V.omissing = true;
2696   do {
2697     Readcard(&cardinfo, &V);
2698     cardhist[(long)cardinfo.art - (long)headercard]++;
2699     switch (cardinfo.art) {
2700
2701     case headercard:
2702     case compndcard:
2703     case sourcecard:
2704     case authorcard:
2705       if (cardhist[(long)cardinfo.art - (long)headercard] == 1) {
2706         for (i = 0; i <= 127; i++)
2707           putc(cardinfo.UU.z[i], tapeout);
2708         putc('\n', tapeout);
2709       }
2710       break;
2711
2712     case ssbondcard:
2713       nss++;
2714       for (i = 0; i <= 1; i++)
2715         StrCopy(ssbonds[nss - 1][i], cardinfo.UU.r[i], sizeof(char6));
2716       break;
2717
2718     case atomcard:
2719       if (cardinfo.UU.U5.residuename != '-' &&
2720           (cardinfo.UU.U5.altloc == 'A' || cardinfo.UU.U5.altloc == ' ')) {
2721         if (strncmp(V.resinfo.aaident, cardinfo.UU.U5.reseqnum, sizeof(char6))) {
2722           Putresidue(&V);
2723           V.resinfo.atompointer = V.latom;
2724           StrCopy(V.resinfo.aaident, cardinfo.UU.U5.reseqnum, sizeof(char6));
2725           V.resinfo.aa = cardinfo.UU.U5.residuename;
2726           StrCopy(V.resinfo.threelettercode, cardinfo.UU.U5.aaname,
2727                  sizeof(char4));
2728         }
2729         Getresidue(cardinfo.UU.U5.atomname, cardinfo.UU.U5.coordinates, &V);
2730       }
2731       if ((cardinfo.UU.U5.residuename == '-') && bVerbose) {
2732         printf(" !!! RESIDUE ");
2733         for (i = 0; i <= 3; i++)
2734           putchar(cardinfo.UU.U5.aaname[i]);
2735         for (i = 0; i <= 5; i++)
2736           putchar(cardinfo.UU.U5.reseqnum[i]);
2737         printf(" HAS NONSTANDARD NAME.\n");
2738         printf("     RESIDUE WILL BE ");
2739         printf("IGNORED !!!\n");
2740       }
2741       if ((cardinfo.UU.U5.altloc != 'A' && cardinfo.UU.U5.altloc != ' ') &&
2742           bVerbose) {
2743         printf(" !!! IN RESIDUE");
2744         for (i = 0; i <= 3; i++)
2745           printf(" %c", cardinfo.UU.U5.aaname[i]);
2746         for (i = 0; i <= 5; i++)
2747           putchar(cardinfo.UU.U5.reseqnum[i]);
2748         printf(" ALTERNATE LOCATION INDICATOR ");
2749         printf("IS %c AND\n", cardinfo.UU.U5.altloc);
2750         printf("     NOT BLANK OR A. ATOM ");
2751         printf("NAMED ");
2752         for (i = 0; i <= 3; i++)
2753           putchar(cardinfo.UU.U5.atomname[i]);
2754         printf(" WILL BE IGNORED !!!\n\n");
2755       }
2756       break;
2757
2758     case tercard:
2759       Putresidue(&V);
2760       break;
2761
2762     case endcard:
2763       finish = true;
2764       Putresidue(&V);
2765       break;
2766
2767     case othercard:
2768       /* blank case */
2769       break;
2770     }
2771   } while (!(V.corelimit || finish));
2772   if ((V.corelimit) && bVerbose) {
2773     printf(" !!! NUMBER OF ATOMS OR RESIDUES EXCEEDS ");
2774     printf("STORAGE CAPACITY !!!\n");
2775   }
2776   if (!Nochainbreak(*V.lchain, *V.lchain))
2777     (*V.lchain)--;
2778   if ((V.hatoms > 0) && bVerbose) {
2779     printf(" !!! %12ld HYDROGEN OR DEUTERIUM ATOMS WERE IGNORED\n", V.hatoms);
2780     printf("     IN THE CALCULATION OF SIDE CHAIN SOLVENT \n");
2781     printf("     ACCESSIBILITY !!!\n");
2782   }
2783   if (bVerbose) {
2784     if (cardhist[0] < 1)
2785       printf(" !!! HEADER-CARD MISSING !!!\n");
2786     if (cardhist[(long)compndcard - (long)headercard] < 1)
2787       printf(" !!! COMPOUND-CARD MISSING !!!\n");
2788     if (cardhist[(long)sourcecard - (long)headercard] < 1)
2789       printf(" !!! SOURCE-CARD MISSING !!!\n");
2790     if (cardhist[(long)authorcard - (long)headercard] < 1)
2791       printf(" !!! AUTHOR CARD MISSING !!!\n");
2792     if (*V.lchain < 1) {
2793       printf(" !!! NO RESIDUE WITH COMPLETE BACKBONE !!!\n");
2794       longjmp(_JL99, 1);
2795     }
2796     if (V.latom == 0)
2797       printf(" !!! ALL SIDECHAIN COORDINATES MISSING !!!\n");
2798   }
2799 }  /* Inputcoordinates */
2800
2801 #undef MAXSIDEATOMS
2802
2803
2804 /***/
2805 /*--------------------------------------------------------------------*/
2806
2807 Static boolean Testbond(i, j)
2808 long i, j;
2809 {
2810   /* TESTBOND IS TRUE IF I IS DONOR[=NH] TO J, OTHERWISE FALSE */
2811   backbone *WITH;
2812
2813   WITH = &chain[i];
2814   return (WITH->acceptor[0].residue == j && WITH->acceptor[0].energy < HBHIGH ||
2815           WITH->acceptor[1].residue == j && WITH->acceptor[1].energy < HBHIGH);
2816 }  /* Testbond */
2817
2818
2819 /***/
2820
2821 Local boolean Testssbond(i, j)
2822 long i, j;
2823 {
2824   boolean ssbond;
2825   long k;
2826
2827   ssbond = false;
2828   k = 1;
2829   if (!(Nochainbreak(i, i) & Nochainbreak(j, j)))
2830     return ssbond;
2831   while (!(ssbond || k > nss)) {
2832     ssbond = ((!strncmp(chain[i].aaident, ssbonds[k - 1][0], sizeof(char6)) &&
2833                !strncmp(chain[j].aaident, ssbonds[k - 1][1], sizeof(char6))) ||
2834               (!strncmp(chain[i].aaident, ssbonds[k - 1][1], sizeof(char6)) &&
2835                !strncmp(chain[j].aaident, ssbonds[k - 1][0], sizeof(char6))));
2836     k++;
2837   }
2838   return ssbond;
2839 }  /* Testssbond */
2840
2841
2842 /***/
2843 /*--------------------------------------------------------------------*/
2844
2845 Static Void Flagssbonds()
2846 {
2847   boolean ssbond;
2848   Char cc;
2849   long i, j, ii, jj;
2850   double d;
2851   long FORLIM;
2852   backbone *WITH;
2853   long FORLIM1;
2854
2855   /***/
2856
2857   nssintra = 0;
2858   nssinter = 0;
2859   cc = '`';
2860   FORLIM = lchain - 2;
2861   for (i = 1; i <= FORLIM; i++) {
2862     if (chain[i].aa == 'C' && chain[i].nsideatoms > 1) {
2863       ii = chain[i].atompointer + 2;
2864       j = i + 1;
2865       do {
2866         j++;
2867         ssbond = false;
2868         if (chain[j].nsideatoms > 1 && chain[j].aa == 'C')
2869           jj = chain[j].atompointer + 2;
2870         else
2871           jj = 0;
2872         if (jj > 0)
2873           ssbond = (Distance(sidechain[ii - 1], sidechain[jj - 1]) < SSDIST);
2874       } while (!(ssbond || j == lchain));
2875       if (ssbond & (!Testssbond(i, j)))
2876         if (bVerbose) {
2877           printf(" !!! ADDITIONAL SSBOND FOUND BETWEEN ");
2878           printf("RESIDUES ");
2879           Writeresidue(chain[i]);
2880           printf(" AND ");
2881           Writeresidue(chain[j]);
2882           printf(" !!!\n\n");
2883         }
2884     }
2885   }
2886   if (nss > 0) {
2887     FORLIM = lchain - 2;
2888     for (i = 1; i <= FORLIM; i++) {
2889       WITH = &chain[i];
2890       if (WITH->aa == 'C') {
2891         FORLIM1 = lchain;
2892         for (j = i + 2; j <= FORLIM1; j++) {
2893           if (chain[j].aa == 'C') {
2894             if (Testssbond(i, j)) {
2895               if (cc == 'z') {
2896                 if (bVerbose)
2897                   printf(" !!! SS-BRIDGE LABEL RESTART AT a !!!\n");
2898                 cc = '`';
2899               }
2900               cc++;
2901               WITH->aa = cc;
2902               chain[j].aa = cc;
2903               if (Nochainbreak(i, j))
2904                 nssintra++;
2905               else
2906                 nssinter++;
2907               if (WITH->nsideatoms > 1) {
2908                 if (chain[j].nsideatoms > 1) {
2909                   jj = chain[j].atompointer + 2;
2910                   ii = WITH->atompointer + 2;
2911                   d = Distance(sidechain[ii - 1], sidechain[jj - 1]);
2912                   if ((d > SSDIST) && bVerbose) {
2913                     printf(" !!! SSBOND DISTANCE IS%5.1f BETWEEN RESIDUES", d);
2914                     Writeresidue(chain[i]);
2915                     printf(" AND ");
2916                     Writeresidue(chain[j]);
2917                     printf(" !!!\n\n");
2918                   }
2919                 }
2920               }
2921             }
2922           }
2923         }
2924       }
2925     }
2926   }
2927   if ((nss != nssintra + nssinter) && bVerbose)
2928     printf(" !!! ERROR IN SSBOND DATA RECORDS !!!\n");
2929 }  /* Flagssbonds */
2930
2931
2932 /***/
2933 /*--------------------------------------------------------------------*/
2934
2935 Static Void Flagchirality()
2936 {
2937   long i;
2938   double ckap, skap;
2939   long FORLIM;
2940   backbone *WITH;
2941
2942   FORLIM = lchain - 2;
2943   for (i = 2; i <= FORLIM; i++) {
2944     WITH = &chain[i];
2945     if (Nochainbreak(i - 1, i + 2)) {
2946       WITH->alpha = Dihedralangle(chain[i - 1].ca, WITH->ca, chain[i + 1].ca,
2947                                   chain[i + 2].ca);
2948       if (WITH->alpha < 0.0)
2949         WITH->ss[(long)chirality - (long)symbol] = '-';
2950       else
2951         WITH->ss[(long)chirality - (long)symbol] = '+';
2952     }
2953   }
2954   FORLIM = lchain - 2;
2955   /***/
2956   for (i = 3; i <= FORLIM; i++) {
2957     WITH = &chain[i];
2958     if (Nochainbreak(i - 2, i + 2)) {
2959       ckap = Cosangle(chain[i].ca, chain[i - 2].ca, chain[i + 2].ca,
2960                       chain[i].ca);
2961       skap = sqrt(1.0 - ckap * ckap);
2962       WITH->kappa = RADIAN * Atan2(skap, ckap);
2963     }
2964   }
2965 }  /* Flagchirality */
2966
2967
2968 /***/
2969
2970 Local long Bondenergy(i, j)
2971 long i, j;
2972 {
2973   /*RESIDUE I IS DONOR[=NH],J IS ACCEPTOR[=CO] OF THE PROTON IN THE
2974      HYDROGEN BOND. THE BONDENERGY IS IN CAL/MOL */
2975   double dho, dhc, dnc, dno;
2976   long hbe;
2977   backbone *WITH;
2978
2979   hbe = 0;
2980   WITH = &chain[i];
2981   if (WITH->aa == 'P')
2982     return hbe;
2983   dho = Distance(WITH->h, chain[j].o);
2984   dhc = Distance(WITH->h, chain[j].c);
2985   dnc = Distance(WITH->n, chain[j].c);
2986   dno = Distance(WITH->n, chain[j].o);
2987   if (dho < DIST || dhc < DIST || dnc < DIST || dno < DIST)
2988     hbe = HBLOW;
2989   else
2990     hbe = (long)floor(Q / dho - Q / dhc + Q / dnc - Q / dno + 0.5);
2991   if (hbe > HBLOW)
2992     return hbe;
2993   if (bVerbose) {
2994     printf(" !!! CONTACT BETWEEN RESIDUES ");
2995     Writeresidue(chain[i]);
2996     printf(" AND ");
2997     Writeresidue(chain[j]);
2998     printf("  TOO CLOSE !!!\n");
2999   }
3000   hbe = HBLOW;
3001   return hbe;
3002 }  /* Bondenergy */
3003
3004 /***/
3005
3006 Local Void Updatebonds(b, hb)
3007 hydrogenbond *b;
3008 hydrogenbond hb;
3009 {
3010   if (hb.energy < b[0].energy) {
3011     b[1] = b[0];
3012     b[0] = hb;
3013   } else if (hb.energy < b[1].energy)
3014     b[1] = hb;
3015 }  /* Updatebonds */
3016
3017 /***/
3018
3019 Local Void Setbonds(i, j)
3020 long i, j;
3021 {
3022   /*I IS NH, J IS CO*/
3023   hydrogenbond hb;
3024
3025   hb.energy = Bondenergy(i, j);
3026   hb.residue = j;
3027   /* CO(J) IS ACCEPTOR OF NH(I) */
3028   Updatebonds(chain[i].acceptor, hb);
3029   hb.residue = i;
3030   Updatebonds(chain[j].donor, hb);
3031 }  /* Setbond */
3032
3033
3034 /***/
3035 /*--------------------------------------------------------------------*/
3036
3037 Static Void Flaghydrogenbonds()
3038 {
3039   long i, j, FORLIM;
3040   backbone *WITH;
3041   long FORLIM1;
3042
3043   /***/
3044
3045   FORLIM = lchain;
3046   for (i = 1; i <= FORLIM; i++) {
3047     if (Nochainbreak(i, i)) {
3048       WITH = &chain[i];
3049       FORLIM1 = lchain;
3050       for (j = i + 1; j <= FORLIM1; j++) {
3051         if (Nochainbreak(j, j)) {
3052           if (Distance(WITH->ca, chain[j].ca) < CADIST) {
3053             Setbonds(i, j);
3054             if (j != i + 1)
3055               Setbonds(j, i);
3056           }
3057         }
3058       }
3059     }
3060   }
3061 }  /* Flaghydrogenbonds */
3062
3063
3064 /***/
3065
3066 Local Void Ladder(i, j, b)
3067 long i, j;
3068 bridgetyp b;
3069 {
3070   long k;
3071   boolean found;
3072   bridge *WITH;
3073
3074   found = false;
3075   k = 1;
3076   if (b == nobridge || i >= j)
3077     return;
3078   do {
3079     WITH = &bridgetable[k - 1];
3080     if (WITH->ib == 0) {
3081       WITH->ib = i;
3082       WITH->ie = i;
3083       WITH->jb = j;
3084       WITH->je = j;
3085       WITH->from = 0;
3086       WITH->towards = 0;
3087       WITH->btyp = b;
3088       nbridge++;
3089       found = true;
3090     } else {
3091       found = (WITH->btyp == b && i == WITH->ie + 1) & Nochainbreak(WITH->ie,
3092                 i) & (((j == WITH->je + 1 && b == parallel) &
3093                        Nochainbreak(WITH->je, j)) | ((j == WITH->jb - 1 &&
3094                           b == antiparallel) & Nochainbreak(j, WITH->jb)));
3095 /* p2c: dssp.p, line 1609: Note:
3096  * Line breaker spent 1.1+0.26 seconds, 3126 tries on line 1540 [251] */
3097       if (found) {
3098         WITH->ie++;
3099         if (b == parallel)
3100           WITH->je++;
3101         else
3102           WITH->jb--;
3103       } else {
3104         k++;
3105         if (k > MAXBRIDGE) {
3106           if (bVerbose)
3107             printf(" !!! BRIDGETABLE OVERFLOW !!!\n");
3108           longjmp(_JL99, 1);
3109         }
3110       }
3111     }
3112   } while (!found);   /* Ladder */
3113 }
3114
3115 /***/
3116
3117 Local Void Testbridge(i)
3118 long i;
3119 {
3120   long j1, j2, j;
3121   bridgetyp b;
3122
3123   /***/
3124
3125   j1 = 0;
3126   j2 = 0;
3127   j = i + 3;
3128   if (!Nochainbreak(i - 1, i + 1))
3129     return;
3130   while (j2 == 0 && j < lchain) {
3131     if (Nochainbreak(j - 1, j + 1)) {
3132       if ((Testbond(i + 1, j) & Testbond(j, i - 1)) |
3133           (Testbond(j + 1, i) & Testbond(i, j - 1)))
3134         b = parallel;
3135       else if ((Testbond(i + 1, j - 1) & Testbond(j + 1, i - 1)) |
3136                (Testbond(j, i) & Testbond(i, j)))
3137         b = antiparallel;
3138       else
3139         b = nobridge;
3140       if (b != nobridge) {
3141         if (j1 == 0) {
3142           j1 = j;
3143           Ladder(i, j, b);
3144         } else if (j != j1) {
3145           j2 = j;
3146           Ladder(i, j, b);
3147         }
3148       }
3149     }
3150     j++;
3151   }
3152 }  /* Testbridge */
3153
3154 /***/
3155
3156 Local Void Extendladder()
3157 {
3158   long i, j, ib1, jb1, je1;
3159   boolean bulge;
3160   long FORLIM;
3161   bridge *WITH;
3162   long SET[11];
3163
3164   FORLIM = nbridge;
3165   for (i = 1; i <= FORLIM; i++) {
3166     WITH = &bridgetable[i - 1];
3167     j = i + 1;
3168     while (j <= nbridge && WITH->towards == 0) {
3169       ib1 = bridgetable[j - 1].ib;
3170       jb1 = bridgetable[j - 1].jb;
3171       je1 = bridgetable[j - 1].je;
3172       bulge = (Nochainbreak(WITH->ie, ib1) && ib1 - WITH->ie < 6 &&
3173                bridgetable[j - 1].btyp == WITH->btyp &&
3174                bridgetable[j - 1].from == 0);
3175       if (bulge) {
3176         switch (WITH->btyp) {
3177
3178         case parallel:
3179           bulge = (jb1 - WITH->je < 6 && ib1 - WITH->ie < 3 ||
3180                    jb1 - WITH->je < 3) & Nochainbreak(WITH->je, jb1);
3181           break;
3182
3183         case antiparallel:
3184           bulge = (WITH->jb - je1 < 6 && ib1 - WITH->ie < 3 ||
3185                    WITH->jb - je1 < 3) & Nochainbreak(je1, WITH->jb);
3186           break;
3187         }
3188       }
3189       if (bulge) {
3190         WITH->towards = j;
3191         bridgetable[j - 1].from = i;
3192       }
3193       j++;
3194     }
3195   }
3196   FORLIM = nbridge;
3197   for (i = 1; i <= FORLIM; i++) {
3198     WITH = &bridgetable[i - 1];
3199     if (WITH->from == 0) {
3200       P_expset(WITH->linkset, 0L);
3201       j = i;
3202       do {
3203         P_addset(WITH->linkset, (int)j);
3204         j = bridgetable[j - 1].towards;
3205       } while (j != 0);
3206       j = WITH->towards;
3207       while (j != 0) {
3208         P_setcpy(bridgetable[j - 1].linkset, WITH->linkset);
3209         j = bridgetable[j - 1].towards;
3210       }
3211     }
3212   }
3213 }  /* Extendladder */
3214
3215 /* Local variables for Sheet: */
3216 struct LOC_Sheet {
3217   long ladderset[MAXBRIDGE / 32 + 2], sheetset[MAXBRIDGE / 32 + 2];
3218 } ;
3219
3220 /***/
3221
3222 Local boolean Link(l1, l2)
3223 long l1, l2;
3224 {
3225   /* LINK IS TRUE IF THERE IS A COMMON RESIDUE IN LADDERS L1 AND L2 */
3226   long ib1, ie1, jb1, je1, ib2, ie2, jb2, je2;
3227
3228   ib1 = bridgetable[l1 - 1].ib;
3229   ie1 = bridgetable[l1 - 1].ie;
3230   jb1 = bridgetable[l1 - 1].jb;
3231   je1 = bridgetable[l1 - 1].je;
3232   ib2 = bridgetable[l2 - 1].ib;
3233   ie2 = bridgetable[l2 - 1].ie;
3234   jb2 = bridgetable[l2 - 1].jb;
3235   je2 = bridgetable[l2 - 1].je;
3236   return (ie1 >= ib2 && ib1 <= ie2 || ie1 >= jb2 && ib1 <= je2 ||
3237           je1 >= ib2 && jb1 <= ie2 || je1 >= jb2 && jb1 <= je2);
3238 }  /* Link */
3239
3240 /***/
3241
3242 Local Void Findsheet(LINK)
3243 struct LOC_Sheet *LINK;
3244 {
3245   long l1, l2;
3246   boolean finish;
3247   long FORLIM, FORLIM1;
3248
3249   /***/
3250
3251   P_expset(LINK->sheetset, 0L);
3252   l1 = 0;
3253   if (*LINK->ladderset != 0L) {
3254     do {
3255       l1++;
3256     } while (!P_inset((int)l1, LINK->ladderset));
3257   }
3258   if (l1 > 0)
3259     P_setcpy(LINK->sheetset, bridgetable[l1 - 1].linkset);
3260   if (l1 <= 0)
3261     return;
3262   do {
3263     finish = true;
3264     FORLIM = nbridge;
3265     for (l1 = 1; l1 <= FORLIM; l1++) {
3266       if (P_inset((int)l1, LINK->sheetset)) {
3267         FORLIM1 = nbridge;
3268         for (l2 = 1; l2 <= FORLIM1; l2++) {
3269           if (P_inset((int)l2, LINK->ladderset)) {
3270             if (Link(l1, l2)) {
3271               P_setunion(LINK->sheetset, LINK->sheetset,
3272                          bridgetable[l2 - 1].linkset);
3273               P_setdiff(LINK->ladderset, LINK->ladderset,
3274                         bridgetable[l2 - 1].linkset);
3275               finish = false;
3276             }
3277           }
3278         }
3279       }
3280     }
3281   } while (!finish);   /* Findsheet */
3282 }
3283
3284 /***/
3285
3286 Local Void Sheet()
3287 {
3288   struct LOC_Sheet V;
3289   long asci, i, j;
3290   Char ccs;
3291   long SET[11];
3292   long FORLIM;
3293   bridge *WITH;
3294
3295   /***/
3296
3297   P_expset(V.ladderset, 0L);
3298   FORLIM = nbridge;
3299   for (i = 1; i <= FORLIM; i++)
3300     P_addset(V.ladderset, (int)i);
3301   ccs = '@';
3302   asci = 64;
3303   while (*V.ladderset != 0L) {
3304     ccs++;
3305     if (ccs > 'z') {
3306       if (bVerbose)
3307         printf(" !!! SHEET LABEL RESTART AT A !!!\n");
3308       ccs = 'A';
3309     }
3310     Findsheet(&V);
3311     FORLIM = nbridge;
3312     for (i = 1; i <= FORLIM; i++) {
3313       WITH = &bridgetable[i - 1];
3314       if (P_inset((int)i, V.sheetset) && WITH->from == 0) {
3315         if (asci == 90) {
3316           if (bVerbose)
3317             printf(" !!! STRAND LABEL RESTART AT A !!!\n");
3318           asci = 64;
3319         }
3320         asci++;
3321         if (WITH->btyp == parallel)
3322           WITH->laddername = (Char)(asci + 32);
3323         else
3324           WITH->laddername = (Char)asci;
3325         WITH->sheetname = ccs;
3326         P_setcpy(WITH->linkset, V.sheetset);
3327         j = WITH->towards;
3328         while (j != 0) {
3329           bridgetable[j - 1].laddername = WITH->laddername;
3330           bridgetable[j - 1].sheetname = WITH->sheetname;
3331           P_setcpy(bridgetable[j - 1].linkset, V.sheetset);
3332           j = bridgetable[j - 1].towards;
3333         }
3334       }
3335     }
3336   }
3337 }  /* Sheet */
3338
3339 /***/
3340
3341 Local Void Markstrands()
3342 {
3343   long i, j, l, ib0, ie0, jb0, je0;
3344   structure beta, betai, betaj;
3345   long iset[(long)beta2 - (long)beta1 + 1][9],
3346        jset[(long)beta2 - (long)beta1 + 1][9];
3347   Char cc;
3348   long FORLIM, FORLIM1;
3349   long SET[9];
3350   bridge *WITH;
3351   backbone *WITH1;
3352   long SET1[9];
3353   long SET2[3];
3354   long SET3[3];
3355
3356   FORLIM = nbridge;
3357   for (i = 1; i <= FORLIM; i++) {
3358     if (bridgetable[i - 1].from == 0) {
3359       j = i;
3360       for (beta = beta1;
3361            (long)beta <= (long)beta2;
3362            beta = (structure)((long)beta + 1)) {
3363         P_setcpy(iset[(long)beta - (long)beta1], P_expset(SET, 0L));
3364         P_setcpy(jset[(long)beta - (long)beta1], P_expset(SET, 0L));
3365       }
3366       ib0 = lchain;
3367       ie0 = 0;
3368       jb0 = lchain;
3369       je0 = 0;
3370       do {
3371         WITH = &bridgetable[j - 1];
3372         FORLIM1 = WITH->ie;
3373         for (l = WITH->ib; l <= FORLIM1; l++) {
3374           WITH1 = &chain[l];
3375           for (beta = beta1;
3376                (long)beta <= (long)beta2;
3377                beta = (structure)((long)beta + 1))
3378             P_setcpy(iset[(long)beta - (long)beta1], P_setunion(SET1,
3379                        iset[(long)beta - (long)beta1],
3380                        P_addset(P_expset(SET, 0L),
3381                                 WITH1->ss[(long)beta - (long)symbol])));
3382         }
3383         FORLIM1 = WITH->je;
3384         for (l = WITH->jb; l <= FORLIM1; l++) {
3385           WITH1 = &chain[l];
3386           for (beta = beta1;
3387                (long)beta <= (long)beta2;
3388                beta = (structure)((long)beta + 1))
3389             P_setcpy(jset[(long)beta - (long)beta1], P_setunion(SET1,
3390                        jset[(long)beta - (long)beta1],
3391                        P_addset(P_expset(SET, 0L),
3392                                 WITH1->ss[(long)beta - (long)symbol])));
3393         }
3394         if (WITH->ib < ib0)
3395           ib0 = WITH->ib;
3396         if (WITH->ie > ie0)
3397           ie0 = WITH->ie;
3398         if (WITH->jb < jb0)
3399           jb0 = WITH->jb;
3400         if (WITH->je > je0)
3401           je0 = WITH->je;
3402         j = WITH->towards;
3403       } while (j != 0);
3404       j = i;
3405       if (P_setequal(iset[0], P_addset(P_expset(SET2, 0L), ' ')))
3406         betai = beta1;
3407       else
3408         betai = beta2;
3409       if (P_setequal(jset[0], P_addset(P_expset(SET2, 0L), ' ')))
3410         betaj = beta1;
3411       else
3412         betaj = beta2;
3413       if ((!P_setequal(iset[(long)betai - (long)beta1],
3414                        P_addset(P_expset(SET2, 0L), ' '))) |
3415           (!P_setequal(jset[(long)betaj - (long)beta1],
3416                        P_addset(P_expset(SET3, 0L), ' '))))
3417         if (bVerbose)
3418           printf(" !!! STRAND COLUMN OVERWRITTEN !!!\n");
3419       do {
3420         WITH = &bridgetable[j - 1];
3421         FORLIM1 = WITH->ie;
3422         for (l = WITH->ib; l <= FORLIM1; l++) {
3423           WITH1 = &chain[l];
3424           WITH1->ss[(long)betai - (long)symbol] = WITH->laddername;
3425           if (WITH->btyp == parallel)
3426             WITH1->partner[(long)betai - (long)beta1] = WITH->jb + l - WITH->ib;
3427           else
3428             WITH1->partner[(long)betai - (long)beta1] = WITH->je - l + WITH->ib;
3429         }
3430         FORLIM1 = WITH->je;
3431         for (l = WITH->jb; l <= FORLIM1; l++) {
3432           WITH1 = &chain[l];
3433           WITH1->ss[(long)betaj - (long)symbol] = WITH->laddername;
3434           if (WITH->btyp == parallel)
3435             WITH1->partner[(long)betaj - (long)beta1] = WITH->ib + l - WITH->jb;
3436           else
3437             WITH1->partner[(long)betaj - (long)beta1] = WITH->ie - l + WITH->jb;
3438         }
3439         j = WITH->towards;
3440       } while (j != 0);
3441       if (ib0 == ie0)
3442         cc = 'B';
3443       else
3444         cc = 'E';
3445       for (j = ib0; j <= ie0; j++) {
3446         WITH1 = &chain[j];
3447         if (WITH1->ss[0] != 'E')
3448           WITH1->ss[0] = cc;
3449       }
3450       for (j = jb0; j <= je0; j++) {
3451         WITH1 = &chain[j];
3452         if (WITH1->ss[0] != 'E')
3453           WITH1->ss[0] = cc;
3454       }
3455     }
3456   }
3457   FORLIM = nbridge;
3458   for (j = 0; j < FORLIM; j++) {
3459     WITH = &bridgetable[j];
3460     FORLIM1 = WITH->ie;
3461     for (l = WITH->ib; l <= FORLIM1; l++)
3462       chain[l].sheetlabel = WITH->sheetname;
3463     FORLIM1 = WITH->je;
3464     for (l = WITH->jb; l <= FORLIM1; l++)
3465       chain[l].sheetlabel = WITH->sheetname;
3466   }
3467 }  /* Markstrands */
3468
3469
3470 /***/
3471 /*--------------------------------------------------------------------*/
3472
3473 Static Void Flagbridge()
3474 {
3475   long i, FORLIM;
3476   bridge *WITH;
3477
3478   /***/
3479
3480   for (i = 0; i < MAXBRIDGE; i++) {
3481     WITH = &bridgetable[i];
3482     WITH->ib = 0;
3483     WITH->ie = 0;
3484     WITH->jb = 0;
3485     WITH->je = 0;
3486     WITH->btyp = nobridge;
3487   }
3488   nbridge = 0;
3489   FORLIM = lchain;
3490   for (i = 2; i < FORLIM; i++)
3491     Testbridge(i);
3492   if (nbridge <= 0)
3493     return;
3494   Extendladder();
3495   Sheet();
3496   Markstrands();
3497 }  /* Flagbridge */
3498
3499
3500 /***/
3501
3502 Local Void Flagsymbol()
3503 {
3504   /* FLAGS ALPHA HELICES AND TURNS IN SYMBOL COLUMN */
3505   long i, j, k;
3506   Char cc;
3507   long nhset[9];
3508   structure turn;
3509   boolean empty;
3510   long FORLIM;
3511   backbone *WITH;
3512
3513   P_addset(P_expset(nhset, 0L), '>');
3514   P_addset(nhset, 'X');
3515   FORLIM = lchain - 4;
3516   for (i = 2; i <= FORLIM; i++) {
3517     if (P_inset(chain[i - 1].ss[(long)turn4 - (long)symbol], nhset) &
3518         P_inset(chain[i].ss[(long)turn4 - (long)symbol], nhset)) {
3519       for (j = i; j <= i + 3; j++)
3520         chain[j].ss[0] = 'H';
3521     }
3522   }
3523   FORLIM = lchain - 3;
3524   for (i = 2; i <= FORLIM; i++) {
3525     if (P_inset(chain[i - 1].ss[(long)turn3 - (long)symbol], nhset) &
3526         P_inset(chain[i].ss[(long)turn3 - (long)symbol], nhset)) {
3527       empty = true;
3528       for (j = i; j <= i + 2; j++) {
3529         WITH = &chain[j];
3530         if (WITH->ss[0] != 'G' && WITH->ss[0] != ' ')
3531           empty = false;
3532       }
3533       if (empty) {
3534         for (j = i; j <= i + 2; j++)
3535           chain[j].ss[0] = 'G';
3536       }
3537     }
3538   }
3539   FORLIM = lchain - 5;
3540   for (i = 2; i <= FORLIM; i++) {
3541     if (P_inset(chain[i - 1].ss[(long)turn5 - (long)symbol], nhset) &
3542         P_inset(chain[i].ss[(long)turn5 - (long)symbol], nhset)) {
3543       empty = true;
3544       for (j = i; j <= i + 4; j++) {
3545         WITH = &chain[j];
3546         if (WITH->ss[0] != 'I' && WITH->ss[0] != ' ')
3547           empty = false;
3548       }
3549       if (empty) {
3550         for (j = i; j <= i + 4; j++)
3551           chain[j].ss[0] = 'I';
3552       }
3553     }
3554   }
3555   FORLIM = lchain;
3556   for (i = 2; i < FORLIM; i++) {
3557     WITH = &chain[i];
3558     if (WITH->ss[0] == ' ') {
3559       cc = ' ';
3560       j = 1;
3561       for (turn = turn3;
3562            (long)turn <= (long)turn5;
3563            turn = (structure)((long)turn + 1)) {
3564         j++;
3565         for (k = 1; k <= j; k++) {
3566           if (i > k) {
3567             if (P_inset(chain[i - k].ss[(long)turn - (long)symbol], nhset))
3568               cc = 'T';
3569           }
3570         }
3571       }
3572       if (cc == ' ')
3573         cc = WITH->ss[(long)bend - (long)symbol];
3574       WITH->ss[0] = cc;
3575     }
3576   }
3577 }  /* Flagsymbol */
3578
3579
3580 /***/
3581 /*--------------------------------------------------------------------*/
3582
3583 Static Void Flagturn()
3584 {
3585   long i, j, k;
3586   structure turn;
3587   Char cc;
3588   long FORLIM1;
3589   backbone *WITH;
3590
3591   /***/
3592
3593   k = 2;
3594   cc = '2';
3595   for (turn = turn3; (long)turn <= (long)turn5; turn = (structure)((long)turn + 1)) {
3596     k++;
3597     cc++;
3598     FORLIM1 = lchain - k;
3599     for (i = 1; i <= FORLIM1; i++) {
3600       if (Nochainbreak(i, i + k)) {
3601         if (Testbond(i + k, i)) {
3602           chain[i + k].ss[(long)turn - (long)symbol] = '<';
3603           for (j = 1; j < k; j++) {
3604             WITH = &chain[i + j];
3605             if (WITH->ss[(long)turn - (long)symbol] == ' ')
3606               WITH->ss[(long)turn - (long)symbol] = cc;
3607           }
3608           WITH = &chain[i];
3609           if (WITH->ss[(long)turn - (long)symbol] == '<')
3610             WITH->ss[(long)turn - (long)symbol] = 'X';
3611           else
3612             WITH->ss[(long)turn - (long)symbol] = '>';
3613         }
3614       }
3615     }
3616   }
3617   FORLIM1 = lchain;
3618   for (i = 1; i <= FORLIM1; i++) {
3619     WITH = &chain[i];
3620     if (WITH->kappa != 360.0 && WITH->kappa > 70.0)
3621       WITH->ss[(long)bend - (long)symbol] = 'S';
3622   }
3623   Flagsymbol();
3624 }  /* Flagturn */
3625
3626
3627 /* Local variables for Flagaccess: */
3628 struct LOC_Flagaccess {
3629   long np;
3630   vector p[NFACE];
3631   double wp[NFACE];
3632 } ;
3633
3634 /* Local variables for Polyeder: */
3635 struct LOC_Polyeder {
3636   struct LOC_Flagaccess *LINK;
3637 } ;
3638
3639 /***/
3640
3641 Local Void Triangle(x1, x2, x3, level, LINK)
3642 double *x1, *x2, *x3;
3643 long level;
3644 struct LOC_Polyeder *LINK;
3645 {
3646   long k, level1;
3647   double xnorm;
3648   vector x4, x5, x6;
3649
3650   if (level > 0) {
3651     level1 = level - 1;
3652     for (k = 0; k <= 2; k++) {
3653       x4[k] = x1[k] + x2[k];
3654       x5[k] = x2[k] + x3[k];
3655       x6[k] = x1[k] + x3[k];
3656     }
3657     Norm(x4, &xnorm);
3658     Norm(x5, &xnorm);
3659     Norm(x6, &xnorm);
3660     Triangle(x1, x4, x6, level1, LINK);
3661     Triangle(x4, x2, x5, level1, LINK);
3662     Triangle(x4, x5, x6, level1, LINK);
3663     Triangle(x5, x3, x6, level1, LINK);
3664     return;
3665   }
3666   for (k = 0; k <= 2; k++)
3667     x6[k] = x1[k] + x2[k] + x3[k];
3668   Norm(x6, &xnorm);
3669   LINK->LINK->np++;
3670   VecCopy(LINK->LINK->p[LINK->LINK->np - 1], x6);
3671   Diff(x3, x1, x5);
3672   Diff(x2, x1, x4);
3673   Cross(x5, x4, x6);
3674   Norm(x6, &xnorm);
3675   LINK->LINK->wp[LINK->LINK->np - 1] = xnorm / 2.0;
3676 }  /* Triangle */
3677
3678 /***/
3679
3680 Local Void Polyeder(LINK)
3681 struct LOC_Flagaccess *LINK;
3682 {  /* GENERATES ALL 12 VERTICES OF ICOSAHEDRON */
3683   struct LOC_Polyeder V;
3684   vector v[12];
3685   double a, b;
3686   long i, j, k, level, FORLIM;
3687
3688   /***/
3689
3690   V.LINK = LINK;
3691   k = 0;
3692   a = YVERTEX;
3693   b = ZVERTEX;
3694   for (i = 1; i <= 2; i++) {
3695     a = -a;
3696     for (j = 1; j <= 2; j++) {
3697       b = -b;
3698       k++;
3699       v[k - 1][0] = 0.0;
3700       v[k - 1][1] = a;
3701       v[k - 1][2] = b;
3702       k++;
3703       v[k - 1][0] = b;
3704       v[k - 1][1] = 0.0;
3705       v[k - 1][2] = a;
3706       k++;
3707       v[k - 1][0] = a;
3708       v[k - 1][1] = b;
3709       v[k - 1][2] = 0.0;
3710     }
3711   }
3712   LINK->np = 0;
3713   level = ORDER;
3714   /* GET ALL 20 FACES OF ICOSAHEDRON */
3715   for (i = 0; i <= 9; i++) {   /* FIND INTEGRATION POINTS */
3716     for (j = i + 1; j <= 10; j++) {
3717       if (Distance(v[i], v[j]) < 1.1) {
3718         for (k = j + 1; k <= 11; k++) {
3719           if ((Distance(v[i], v[k]) < 1.1) & (Distance(v[j], v[k]) < 1.1))
3720             Triangle(v[i], v[j], v[k], level, &V);
3721         }
3722       }
3723     }
3724   }
3725   a = 0.0;
3726   FORLIM = LINK->np;
3727   for (i = 0; i < FORLIM; i++)
3728     a += LINK->wp[i];
3729   a = FOURPI / a;
3730   FORLIM = LINK->np;
3731   for (i = 0; i < FORLIM; i++)
3732     LINK->wp[i] *= a;
3733 }  /* Polyeder (enurD idu) */
3734
3735 /* Local variables for Surface: */
3736 struct LOC_Surface {
3737   struct LOC_Flagaccess *LINK;
3738   long nx;
3739   vector x[MAXPACK];
3740   double rx[MAXPACK];
3741 } ;
3742
3743 /***/
3744
3745 Local boolean Step(xx, LINK)
3746 double *xx;
3747 struct LOC_Surface *LINK;
3748 {
3749   long k;
3750   boolean one;
3751   double TEMP;
3752
3753   one = true;
3754   k = 1;
3755   while (k <= LINK->nx && one) {
3756     TEMP = LINK->rx[k - 1] + RWATER;
3757     if (Distsq(xx, LINK->x[k - 1]) < TEMP * TEMP)
3758       one = false;
3759     else
3760       k++;
3761   }
3762   return one;
3763 }  /* Step */
3764
3765 /* Local variables for Liste: */
3766 struct LOC_Liste {
3767   struct LOC_Surface *LINK;
3768 } ;
3769
3770 /***/
3771
3772 Local Void Listentry(xx, yy, d, r, LINK)
3773 double *xx, *yy;
3774 double d, r;
3775 struct LOC_Liste *LINK;
3776 {
3777   vector zz;
3778   double delta;
3779
3780   delta = Distance(xx, yy);
3781   if (delta >= d + r)
3782     return;
3783   if (delta <= EPS)
3784     return;
3785   LINK->LINK->nx++;
3786   if (LINK->LINK->nx > MAXPACK) {
3787     if (bVerbose)
3788       printf(" !!! TABLE OVERFLOW IN FLAGACCESS !!!\n");
3789     longjmp(_JL99, 1);
3790     return;
3791   }
3792   Diff(yy, xx, zz);
3793   VecCopy(LINK->LINK->x[LINK->LINK->nx - 1], zz);
3794   LINK->LINK->rx[LINK->LINK->nx - 1] = r;
3795 }  /* Listentry */
3796
3797 /***/
3798
3799 Local Void Liste(xx, rxx, LINK)
3800 double *xx;
3801 double rxx;
3802 struct LOC_Surface *LINK;
3803 {
3804   struct LOC_Liste V;
3805   long i, k;
3806   double d;
3807   long FORLIM;
3808   backbone *WITH;
3809   long FORLIM1;
3810
3811   /***/
3812
3813   V.LINK = LINK;
3814   LINK->nx = 0;
3815   d = rxx + RWATER + RWATER;
3816   FORLIM = lchain;
3817   for (i = 1; i <= FORLIM; i++) {
3818     if (Nochainbreak(i, i)) {
3819       WITH = &chain[i];
3820       if (Distance(xx, WITH->ca) < d + RESRAD) {
3821         Listentry(xx, WITH->n, d, RN, &V);
3822         Listentry(xx, WITH->ca, d, RCA, &V);
3823         Listentry(xx, WITH->c, d, RC, &V);
3824         Listentry(xx, WITH->o, d, RO, &V);
3825         if (WITH->nsideatoms > 0) {
3826           FORLIM1 = WITH->nsideatoms;
3827           for (k = 0; k < FORLIM1; k++)
3828             Listentry(xx, sidechain[WITH->atompointer + k], d, RSIDEATOM, &V);
3829         }
3830       }
3831     }
3832   }
3833 }  /* Liste */
3834
3835 /***/
3836
3837 Local double Surface(xatom, ratom, LINK)
3838 double *xatom;
3839 double ratom;
3840 struct LOC_Flagaccess *LINK;
3841 {
3842   struct LOC_Surface V;
3843   long i, j;
3844   double f, radius;
3845   vector xx;
3846   long FORLIM;
3847
3848   /***/
3849
3850   V.LINK = LINK;
3851   Liste(xatom, ratom, &V);
3852   radius = ratom + RWATER;
3853   f = 0.0;
3854   FORLIM = LINK->np;
3855   for (i = 0; i < FORLIM; i++) {
3856     for (j = 0; j <= 2; j++)
3857       xx[j] = LINK->p[i][j] * radius;
3858     if (Step(xx, &V))
3859       f += LINK->wp[i];
3860   }
3861   return (radius * radius * f);
3862 }  /* Surface */
3863
3864
3865 /***/
3866 /*--------------------------------------------------------------------*/
3867
3868 Static Void Flagaccess()
3869 {
3870   struct LOC_Flagaccess V;
3871   long i, k;
3872   double f;
3873   long FORLIM;
3874   backbone *WITH;
3875   long FORLIM1;
3876
3877   /***/
3878
3879   Polyeder(&V);
3880   FORLIM = lchain;
3881   for (i = 1; i <= FORLIM; i++) {
3882     if (Nochainbreak(i, i)) {
3883       WITH = &chain[i];
3884       f = Surface(WITH->n, RN, &V) + Surface(WITH->ca, RCA, &V) +
3885           Surface(WITH->c, RC, &V) + Surface(WITH->o, RO, &V);
3886       if (WITH->nsideatoms > 0) {
3887         FORLIM1 = WITH->nsideatoms;
3888         for (k = 0; k < FORLIM1; k++)
3889           f += Surface(sidechain[WITH->atompointer + k], RSIDEATOM, &V);
3890       }
3891       WITH->access = (long)floor(f + 0.5);
3892     }
3893   }
3894 }  /* Flagaccess */
3895
3896
3897 /***/
3898
3899 Local Void Statistics()
3900 {
3901   long i, j, k, nchain, nres, nhbond, lhelix;
3902   bridgetyp b;
3903   Char cc;
3904   double Surface;
3905   long nhbturn[11];
3906   long ladderset[MAXBRIDGE / 32 + 2];
3907   long hbridge[(long)antiparallel - (long)parallel + 1];
3908   long helixhist[MAXHIST], sheethist[MAXHIST];
3909   long betahist[(long)antiparallel - (long)parallel + 1][MAXHIST];
3910   long FORLIM, FORLIM1;
3911   backbone *WITH;
3912   bridge *WITH1;
3913   long SET[11];
3914   long SET1[257];
3915
3916   lhelix = 0;
3917   nhbond = 0;
3918   nchain = 0;
3919   nres = 0;
3920   for (i = 0; i < MAXHIST; i++) {
3921     for (b = parallel; (long)b <= (long)antiparallel; b = (bridgetyp)((long)b + 1))
3922       betahist[(long)b - (long)parallel][i] = 0;
3923     helixhist[i] = 0;
3924     sheethist[i] = 0;
3925   }
3926   Surface = 0.0;
3927   for (k = 0; k <= 10; k++)
3928     nhbturn[k] = 0;
3929   for (b = parallel; (long)b <= (long)antiparallel; b = (bridgetyp)((long)b + 1))
3930     hbridge[(long)b - (long)parallel] = 0;
3931   FORLIM = lchain;
3932   for (i = 0; i <= FORLIM; i++) {
3933     WITH = &chain[i];
3934     if (Nochainbreak(i, i)) {
3935       nres++;
3936       Surface += WITH->access;
3937       for (j = 0; j <= 1; j++) {
3938         if (WITH->donor[j].energy < HBHIGH) {
3939           nhbond++;
3940           k = WITH->donor[j].residue - i;
3941           if (labs(k) < 6)
3942             nhbturn[k + 5]++;
3943         }
3944       }
3945     } else
3946       nchain++;
3947     if (WITH->ss[0] == 'H')
3948       lhelix++;
3949     else if (lhelix > 0) {
3950       if (lhelix > MAXHIST)
3951         lhelix = MAXHIST;
3952       helixhist[lhelix - 1]++;
3953       lhelix = 0;
3954     }
3955   }
3956   if (nbridge > 0) {
3957     FORLIM = nbridge;
3958     for (i = 1; i <= FORLIM; i++) {
3959       WITH1 = &bridgetable[i - 1];
3960       hbridge[(long)WITH1->btyp - (long)parallel] += WITH1->ie - WITH1->ib + 2;
3961       if (WITH1->from == 0) {
3962         j = i;
3963         k = 0;
3964         do {
3965           k += bridgetable[j - 1].ie - bridgetable[j - 1].ib + 1;
3966           j = bridgetable[j - 1].towards;
3967         } while (j != 0);
3968         if (k > MAXHIST)
3969           k = MAXHIST;
3970         betahist[(long)WITH1->btyp - (long)parallel][k - 1]++;
3971       }
3972     }
3973   }
3974   if (nbridge > 0) {
3975     P_expset(ladderset, 0L);
3976     FORLIM = nbridge;
3977     for (i = 1; i <= FORLIM; i++)
3978       P_addset(ladderset, (int)i);
3979     FORLIM = nbridge;
3980     for (i = 1; i <= FORLIM; i++) {
3981       WITH1 = &bridgetable[i - 1];
3982       if ((WITH1->from == 0) & P_inset((int)i, ladderset)) {
3983         if (!P_setequal(P_addset(P_expset(SET1, 0L), (int)i), WITH1->linkset) ||
3984             WITH1->ie > WITH1->ib) {
3985           k = 0;
3986           FORLIM1 = nbridge;
3987           for (j = 1; j <= FORLIM1; j++) {
3988             if ((bridgetable[j - 1].from == 0) & P_inset((int)j, WITH1->linkset))
3989               k++;
3990           }
3991           sheethist[k - 1]++;
3992         }
3993         P_setdiff(ladderset, ladderset, WITH1->linkset);
3994       }
3995     }
3996   }
3997   fprintf(tapeout,
3998     "%5ld%3ld%3ld%3ld%3ld TOTAL NUMBER OF RESIDUES, NUMBER OF CHAINS, NUMBER OF SS-BRIDGES(TOTAL,INTRACHAIN,INTERCHAIN)                .\n",
3999     nres, nchain, nssinter + nssintra, nssintra, nssinter);
4000   fprintf(tapeout,
4001     "%8.1f   ACCESSIBLE SURFACE OF PROTEIN (ANGSTROM**2)                                                                         .\n",
4002     Surface);
4003   fprintf(tapeout,
4004     "%5ld%5.1f   TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(J)  , SAME NUMBER PER 100 RESIDUES                              .\n",
4005     nhbond, 100.0 * nhbond / nres);
4006   i = hbridge[0];
4007   j = hbridge[(long)antiparallel - (long)parallel];
4008   fprintf(tapeout,
4009     "%5ld%5.1f   TOTAL NUMBER OF HYDROGEN BONDS IN     PARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES                              .\n",
4010     i, 100.0 * i / nres);
4011   fprintf(tapeout,
4012     "%5ld%5.1f   TOTAL NUMBER OF HYDROGEN BONDS IN ANTIPARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES                              .\n",
4013     j, 100.0 * j / nres);
4014   for (i = -5; i <= 5; i++) {
4015     if (i < 0)
4016       cc = '-';
4017     else
4018       cc = '+';
4019     k = labs(i);
4020     fprintf(tapeout,
4021       "%5ld%5.1f   TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(I%c%ld), SAME NUMBER PER 100 RESIDUES                              .\n",
4022       nhbturn[i + 5], 100.0 * nhbturn[i + 5] / nres, cc, k);
4023   }
4024   for (i = 1; i <= MAXHIST; i++)
4025     fprintf(tapeout, "%3ld", i);
4026   fprintf(tapeout, "     *** HISTOGRAMS OF ***           .\n");
4027   for (i = 0; i < MAXHIST; i++)
4028     fprintf(tapeout, "%3ld", helixhist[i]);
4029   fprintf(tapeout, "    RESIDUES PER ALPHA HELIX         .\n");
4030   for (i = 0; i < MAXHIST; i++)
4031     fprintf(tapeout, "%3ld", betahist[0][i]);
4032   fprintf(tapeout, "    PARALLEL BRIDGES PER LADDER      .\n");
4033   for (i = 0; i < MAXHIST; i++)
4034     fprintf(tapeout, "%3ld", betahist[(long)antiparallel - (long)parallel][i]);
4035   fprintf(tapeout, "    ANTIPARALLEL BRIDGES PER LADDER  .\n");
4036   for (i = 0; i < MAXHIST; i++)
4037     fprintf(tapeout, "%3ld", sheethist[i]);
4038   fprintf(tapeout, "    LADDERS PER SHEET                .\n");
4039 }  /* Statistics */
4040
4041 /***/
4042
4043 Local Void Writehb(i, hb)
4044 long i;
4045 hydrogenbond hb;
4046 {
4047   double e;
4048
4049   if (hb.residue != 0)
4050     hb.residue -= i;
4051   e = hb.energy / 1000.0;
4052   fprintf(tapeout, "%4ld,%4.1f", hb.residue, e);
4053 }  /* Writehb */
4054
4055
4056 /***/
4057 /*--------------------------------------------------------------------*/
4058
4059 Static Void Printout()
4060 {
4061   long i, j;
4062   structure s;
4063   double phi, psi, tco;
4064   long FORLIM;
4065   backbone *WITH;
4066
4067   /***/
4068
4069   Statistics();
4070   fprintf(tapeout,
4071     "  #  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");
4072   FORLIM = lchain;
4073   for (i = 1; i <= FORLIM; i++) {
4074     WITH = &chain[i];
4075     fprintf(tapeout, "%5ld ", i);
4076     for (j = 0; j <= 5; j++)
4077       putc(WITH->aaident[j], tapeout);
4078     fprintf(tapeout, " %c  %c ", WITH->aa, WITH->ss[0]);
4079     for (s = turn3; (long)s <= (long)beta2; s = (structure)((long)s + 1))
4080       putc(WITH->ss[(long)s - (long)symbol], tapeout);
4081     for (s = beta1; (long)s <= (long)beta2; s = (structure)((long)s + 1))
4082       fprintf(tapeout, "%4ld", WITH->partner[(long)s - (long)beta1]);
4083     fprintf(tapeout, "%c%4ld ", WITH->sheetlabel, WITH->access);
4084     for (j = 0; j <= 1; j++) {
4085       Writehb(i, WITH->acceptor[j]);
4086       Writehb(i, WITH->donor[j]);
4087     }
4088     phi = 360.0;
4089     psi = 360.0;
4090     tco = 0.0;
4091     if (Nochainbreak(i - 1, i)) {
4092       phi = Dihedralangle(chain[i - 1].c, WITH->n, WITH->ca, WITH->c);
4093       tco = Cosangle(WITH->c, WITH->o, chain[i - 1].c, chain[i - 1].o);
4094     }
4095     if (Nochainbreak(i, i + 1))
4096       psi = Dihedralangle(WITH->n, WITH->ca, WITH->c, chain[i + 1].n);
4097     fprintf(tapeout, "%8.3f%6.1f%6.1f%6.1f%6.1f%7.1f%7.1f%7.1f\n",
4098             tco, WITH->kappa, WITH->alpha, phi, psi, WITH->ca[0], WITH->ca[1],
4099             WITH->ca[2]);
4100   }
4101 }  /* Printout */
4102
4103 Static Void Usage()
4104 {
4105   fprintf(stderr,"Usage: dssp [-na] pdb_file dssp_file\n");
4106   fprintf(stderr,"the -na flag disables the calculation of accessible surface\n");
4107 }
4108
4109 /***/
4110 /*--------------------------------------------------------------------*/
4111
4112 void printit()
4113 {
4114   printf(" \n");
4115   printf("                           DSSP\n");
4116   printf("            by Wolfgang Kabsch and Chris Sander\n");
4117   printf("\n");
4118   printf("Defines secondary structure and solvent exposure of proteins from\n");
4119   printf("atomic coordinates as given in Protein Data Bank format.   \n");
4120   printf("\n");
4121   printf("_________________________________________________________________________\n");
4122   printf("This version licensed to Ethan Benatan at Univ_Pittsburgh               \n");
4123   printf("for academic purposes.                                                 \n");
4124   printf("Do not redistribute. \n");
4125   printf("\n");
4126   printf("Commercial licenses available on request.                                \n");
4127   printf("\n");
4128   printf("Copyright by Wolfgang Kabsch and Chris Sander, 1983, 1985, 1988. \n");
4129   printf("Fax: +49-6221-387 306\n");
4130   printf("\n");
4131   printf("Algorithm version October 1988. Refer to Biopolymers 22(1983) 2577-2637.\n");
4132   printf("Do report errors if you find any.\n");
4133   printf("\n");
4134   printf("Email: Sander@embl-heidelberg.de \n");
4135   printf("       Kabsch@embl-heidelberg.de \n");
4136   printf("\n");
4137   printf("_________________________________________________________________________\n");
4138   printf("Related databases and datasets available from the Protein\n");
4139   printf("Design Group at EMBL via anonymous ftp from embl-heidelberg.de:\n");
4140   printf("\n");
4141   printf("pdb_select   Representative set of protein structures.\n");
4142   printf("             By U. Hobohm, C. Sander, M. Scharf and R. Schneider.\n");
4143   printf("             See Protein Science 1, 409-417.\n");
4144   printf("DSSP         Dictionary of secondary structures of proteins. \n");
4145   printf("HSSP         Database of sequence-homology derived protein families.\n");
4146   printf("             By C. Sander and R.Schneider.\n");
4147   printf("             See Proteins 9, 56-68 (1991).\n");
4148   printf("FSSP         Database of protein structure families with \n");
4149   printf("             common folding motifs. \n");
4150   printf("             L.Holm, C. Ouzounis, C. Sander, G.Tuparev, G. Vriend\n");
4151   printf("             See Protein Science, in the press (1992).\n");
4152   printf("In the XSSP databases, there is one dataset for each unique or\n");
4153   printf("             representative PDB protein, e.g., 1PPT.HSSP etc.\n");
4154   printf("\n");
4155   printf("Restrictions:Commercial users must apply for a license. \n");
4156   printf("             Not to be used for classified research.\n");
4157 }
4158
4159 void dssp_main(int bDoAcc, int bSetVerbose)
4160 {
4161   int tt; /*TIMELOCK*/
4162   PASCAL_MAIN(0,NULL);
4163   if (setjmp(_JL99))
4164     goto _L99;
4165
4166   bVerbose=bSetVerbose;
4167
4168   tt=time(0); /*TIMELOCK*/
4169
4170   lchain = 0;
4171   Inputcoordinates(&lchain);
4172   if (!Nochainbreak(1L, lchain))
4173     printf(" !!! POLYPEPTIDE CHAIN INTERRUPTED !!!\n");
4174 /*   printf("INPUTCOORDINATES DONE%12ld\n", lchain); */
4175   Flagssbonds();
4176 /*   printf("FLAGSSBONDS DONE\n"); */
4177   Flagchirality();
4178 /*   printf("FLAGCHIRALITY DONE\n"); */
4179   Flaghydrogenbonds();
4180 /*   printf("FLAGHYDROGENBONDS DONE\n"); */
4181   Flagbridge();
4182 /*   printf("FLAGBRIDGE DONE\n"); */
4183   Flagturn();
4184 /*   printf("FLAGTURN DONE\n"); */
4185   if (bDoAcc) {
4186     Flagaccess();
4187 /*     printf("FLAGACCESS DONE\n"); */
4188   }
4189   Printout();
4190 /*   printf("PRINTOUT DONE\n"); */
4191 _L99:
4192   /*fclose(tapein);
4193     fclose(tapeout);*/
4194   ;
4195 }  /* END OF PROGRAM DSSP */
4196