src/contrib/dsspcore.c

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