C PROTEIN DATA BANK SOURCE CODE LSM C AUTHOR. R.MATELA,R.FLETTERICK C ENTRY DATE. 3/82 UNSUPPORTED C LAST REVISION. 3/82 C PURPOSE. COLOR-CODED ALPHA-CARBON MODELS C LANGUAGE. FORTRAN 77,DEC VAX 11 VMS C WARNING. ALL OCCURRENCES OF $ (EXCEPT IN C WARNING. FORMAT STATEMENTS) MUST BE C WARNING. REPLACED BY UNDERLINE (ASCII C WARNING. 137) BEFORE COMPILATION C C **************************************************************** C * * C * LSM * C * * C * VERSION 2.0 FEB. 25,10982 * C * * C * * C **************************************************************** C C C==>THIS PROGRAM GETS THE ALPHA CARBON ATOMS FOR A GIVEN PROTEIN C==>AND FORCES ALL BOND LENGTHS TO 3.8 ANGSTROMS AT THE EXPENSE C==>OF DIHEDRAL AND BOND ANGLES WHILE MINIMIZING THE SHIFT OF C==>INDIVIDUAL ATOMS. THIS IS A LOCAL PROCESS. C C==>THIS IS A NEW ALGORITHM FOR SHIFTING ALPHA - CARBON ATOMS C==>USING A LAGRANGIAN FUNCTION, JACOBIAN AND NEWTON - RAPHSON C==>ITERATION. C C C==>THIS VERSION WILL HANDLE MULTI-SECTION DATA, I.E. DATA WITH C==>MISSING ATOMS. INDIVIDUAL SECTIONS ARE OPTIMIZED. C C C==>NOTE:THIS VERSION DOES ONE ITERATION(SHIFT) PER ATOM PER CYCLE!!! C C==>NOTE:IN SUBROUTINE CHANGE THERE ARE TWO THINGS TO LOOK OUT FOR: C 1) THE MAXIMUM SHIFT ALLOWED C 2) THE DIRECTION OF THE SHIFT C C==> THE MAGNITUDE OF THE SHIFT IS GIVEN BY THE CONSTRUCTED JACOBIAN. C==> TYPICAL VALUES FOR (MAX$SHIFT) ARE 0.001-->1.5. FOR THE CORRECTION C==> FACTOR (COR$FACT) TYPICAL VALUES ARE 0.01-->0.5. BOTH OF THESE C==> VALUES WILL DEPEND UPON YOUR DATA. C C==> TO TRAP FOR A WRONG DIRECTION (VERY RARE) CHECK IF DETERMINANT C==> IS NEGATIVE. IF IT IS CHANGE THE SIGN OF THE GRADIENT VECTOR C==> AND MULTIPLY THIS VECTOR BY SOME SCALAR FACTOR,E.G. 0.1. C C **************************************************************** C * * C * FILES * C * * C **************************************************************** C C C==>NOTE: IN BROOKHAVEN VERSION COMMON BLOCK FILES LSCOM1.FOR AND LBCOM1.FOR C HAVE BEEN EXPANDED IN THE SOURCE AND ARE NOT NECESSARY. C C C==> INPUT: AP****.DAT (E.G. AP1MBN.DAT) C CONTAINS ALPHA-CARBON X,Y,Z AND RES TYPE DATA C ORIGINAL DATA FROM BROOKHAVEN DATA BANK C C RESCOL.DAT C CONTAINS ONE AND THREE LETTER AMINO ACID CODES C AND SHAPELY COLOR ASSIGNMENTS. CAN BE MODIFIED C FOR SPECIFIC APPLICATIONS C C LSCOM1.FOR C CONTAINS COMMON BLOCK INFORMATION FOR THE SHIFT C PORTION OF THE PROGRAM C C LBCOM1.FOR C CONTAINS COMMON BLOCK INFORMATION FOR THE BEND C PORTION OF THE PROGRAM C C C CHARCOM.FOR C CONTAINS COMMON BLOCK INFORMATION FOR CHARACTER C STRINGS C C==> OUTPUT: AOP****.DAT (E.G. AOP1MBN.DAT) C CONTAINS RELEVANT INFORMATION FROM THE SHIFT C PORTION OF THE PROGRAM. FOR EXAMPLE,THE ORIGINAL C BOND LENGTHS,THE NEW BOND LENGTHS,THE OLD AND NEW C COORDINATES,ETC. C C BOP****.DAT (E.G. BOP1MBN.DAT) C CONTAINS RELEVANT INFORMATION FROM THE BEND C PORTION OF THE PROGRAM. FOR EXAMPLE,THE BEND C AND DIHEDERAL/TORSIONAL ANGLES, NEAREST NEIGHBORS C AND COLOR CODES. THIS IS THE ONLY FILE NEEDED TO C CONSTRUCT A SHAPELY MOLECULAR MODEL C C C **************************************************************** C * * C * SAMPLE COLOR ASSIGNMENTS * C * * C **************************************************************** C C A ALA YELLOW C Q GLN RED C N ASN RED C L LEU YELLOW C G GLY WHITE C K LYS RED C R ARG YELLOW C T THR WHITE C P PRO WHITE C I ILE YELLOW C M MET YELLOW C F PHE YELLOW C Y TYR WHITE C C CYS YELLOW C W TRP WHITE C H HIS RED C E GLU RED C D ASP RED C V VAL YELLOW C S SER WHITE C C C **************************************************************** C * * C * AUTHORS * C * * C **************************************************************** C C C DR. RAYMOND MATELA C THE OPEN UNIVERSITY C FACULTY OF TECHNOLOGY C CENTRE FOR CONFIGURATIONAL STUDIES C MILTON KEYNES,BUCKS MK7 6AA C ENGLAND C C DR. ROBERT FLETTERICK C DEPARTMENT OF BIOCHEMISTRY C UNIVERSITY OF CALIFORNIA C SCHOOL OF MEDICINE C SAN FRANCISCO, CALIFORNIA 94143 C U.S.A C C C C INCLUDE 'CHARCOM.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C INCLUDE 'LSCOM1.FOR' DATA YES/'Y'/ C C **************************************************************** C * * C * VARIABLES * C * * C **************************************************************** C C==> NA THE NUMBER OF ALPHA CARBON ATOMS IN MOLECULE C==> RES 3 LETTER CODES FOR RESPECTIVE RESIDUES C==> AP1 ORIGINAL (GIVEN) X,Y,Z COORDINATES C==> CONV CONVERGENCE VALUE TO TERMINATE RUN C==> NCYC THE MAXIMUM NUMBER OF CYCLES THROUGH MOLECULE C==> TBL THE TARGET BOND LENGTH (I.E. 3.8 ANGSTROMS) C==> DF1 VECTOR OF FIRST DERIVATIVES (GRAD OF FUNCTION) C==> JAC JACOBIAN MATRIX C==> AINV INVERSE MATRIX OF JAC C==> C VECTOR CONTAINING THE COMPUTED CHANGES (SHIFTS)+ LAMBDA'S C==> DET DETERMINANT OF JAC C==> XYZP CURRENT (NEW) X,Y,Z COORDINATES C==> BL1 FOR ATOM I, THE (I,I-1) BOND LENGTH C==> BL2 FOR ATOM I, THE (I,I+1) BOND LENGTH C==> BLEN1 THE NEW INDIVIDUAL (I,I+1) BOND LENGTH C==> NBL A VECTOR CONTAINING THE BLEN1'S FOR THE MOLECULE C==> BLEN THE SUM OF ALL THE NEW BOND LENGTHS C==> SSLDEV THE SUM OF THE SQUARES OF THE NEW BOND LENGTH DEVIATIONS C==> AVGLEN THE AVERAGE NEW BOND LENGTH C==> SSPOS THE SUM OF THE SQUARES OF THE NEW POSITION DEVIATIONS C==> LAM1 LAMBDA 1 C==> LAM2 LAMBDA 2 C==> MAX$SHIFT MAXIMUM SHIFT ALLOWED AS COMPUTED IN SUB. CHANGE C==> COR$FACT THE SCALAR FACTOR BY WHICH C(I),I=1,3 IS MULTIPLIED C C **************************************************************** C * * C * SET UP THE EXCESSIVE SHIFT VALUES * C * * C * NOTE: THESE ARE VERY IMPORTANT VALUES USED IN * C * SUBROUTINE CHANGE. IF TROUBLES ARISE WITH BAD * C * BONDS, TWEEK THESE. * C * * C **************************************************************** C MAX$SHIFT = 1.5D0 COR$FACT = 0.1D0 C C **************************************************************** C * * C * SET UP THE NECESSARY FILES * C * * C **************************************************************** C CALL FILES C C **************************************************************** C * * C * SET UP THE VARIOUS ARRAYS * C * * C **************************************************************** C CALL ARRAY$SETUP C C **************************************************************** C * * C * CHECK IF THE MOLECULE IS IN DISJOINT SECTIONS, THAT * C * IS, THE MOLECULE IS MISSING DATA * C * ALSO CHECH IF NEW RESIDUES HAVE BEEN ADDED. * C * * C **************************************************************** C SECT$NUM = 1 NUM$SHOULD$BE = 0 SECT( SECT$NUM,1 ) = 1 SECT( SECT$NUM,3 ) = ORG$NUM( 1 ) C C==>NOTE: THIS PATCH ASSUMES FIRST A.A. UNIT HAS NUMBER 1. C==> MIGHT HAVE TO RE-WRITE THIS CODE FOR CASES. C DO 11,I=1,NA NUM$SHOULD$BE = ORG$NUM( I-1 )+ 1 IF( I .EQ. 1 ) NUM$SHOULD$BE = 1 C C C **************************************************************** C * * C * * C * TRAP FOR ADDED RESIDUES,E.G. 27, 27A, 27B, ETC. * C * * C **************************************************************** C C C IF( ORG$NUM(I) .EQ. ORG$NUM(I-1)) GOTO 11 IF ( ORG$NUM( I ) .EQ. NUM$SHOULD$BE ) GOTO 11 C C==>NUMBERS DO NOT MATCH - END ONE SECTION,START ANOTHER C C==>END FIRST C SECT ( SECT$NUM,2 ) = I - 1 SECT( SECT$NUM,4 ) = ORG$NUM( I-1 ) TYPE *,(SECT(SECT$NUM,J),J=1,4) C C==>START NEW SECTION C SECT$NUM = SECT$NUM + 1 SECT( SECT$NUM,1 ) = I SECT( SECT$NUM,3 ) = ORG$NUM( I ) C C==>RESET THE NUMBER THAT SHOULD BE TO AGREE WITH ORIG DATA C NUM$SHOULD$BE = SECT( SECT$NUM,3 ) 11 CONTINUE SECT(SECT$NUM,2) = NA SECT(SECT$NUM,4) = ORG$NUM(NA) TYPE *,(SECT(SECT$NUM,J),J=1,4) C C **************************************************************** C * * C * START CYCLE THROUGH MOLECULE AND CHECK DATA FOR * C * GARBAGE BOND LENGTHS. * C * * C **************************************************************** C DO 950,IJK=1,SECT$NUM C C==>CLEAR FLAG FOR SECTIONS WITH LESS THAN 3 ATOMS C LOW$ATOM$FLAG = 0 C C **************************************************************** C * * C * INITIALIZE THE CONVERGENCE VALUE. NOTE THAT THE * C * CONVERGENCE TEST IS ON THE SUM OF THE SQUARES OF * C * THE DEVIATIONS OF BOND LENGTHS. * C * * C **************************************************************** C CONV = 10.D0 TYPE *,IJK,SECT$NUM HIGH$ATOM$SECT = SECT( IJK,2 ) NUM$ATOM$SECT = SECT(IJK,2) - SECT(IJK,1) + 1 LOW$ATOM$SECT = SECT( IJK,1) TYPE *,LOW$ATOM$SECT,HIGH$ATOM$SECT C C C **************************************************************** C * * C * * C * TRAP FOR SECTION WITH LESS THAN THREE ATOMS * C * YES IT HAPPENS, LOOK AT P1TIM !!!! * C * * C **************************************************************** C C==>YES IT HAPPENS,LOOK AT P1TIM !!! C IF( (HIGH$ATOM$SECT - LOW$ATOM$SECT) + 1 .LT. 3 ) THEN TYPE *,'NUMBER OF ATOMS IN SECTION LESS THAN 3 ! ' LOW$ATOM$FLAG = 1 CALL SHIFT$OUT CALL ANGLES CALL ANGLE$OUT GOTO 950 END IF DO 15,I=LOW$ATOM$SECT,HIGH$ATOM$SECT-1 BL1 = DSQRT (( AP1(I,1) - AP1(I+1,1 ))**2 1 + ( AP1(I,2) - AP1(I+1,2 ))**2 1 + ( AP1(I,3) - AP1(I+1,3 ))**2) C C **************************************************************** C * * C * SAVE THE ORIGINAL BOND LENGTHS * C * * C **************************************************************** C ORBL(I) = BL1 IF (BL1 .LT. 2.3D0 .OR. BL1 .GT. 5.3D0) THEN TYPE 610,I,BL1,RES(I),(AP1(I,J),J=1,3) 610 FORMAT(1X,' BAD BOND= ',I4,' LENGTH= ',D15.7,' RES= ', 1 A3,3X,' XYZ =',3(1X,D15.7)) END IF 15 CONTINUE C C **************************************************************** C * * C * SET THE MAXIMUM NUMBER OF CYCLES THROUGH THE * C * MOLECULE. * C * * C **************************************************************** C C NCYC = 10 C C **************************************************************** C * * C * SET THE TARGET (FIXED) BOND LENGTH. * C * * C **************************************************************** C TBL=3.8D0 C C **************************************************************** C * * C * INITIALIZE X, Y, Z PRIME COORDINATES AND CLEAR THE * C * LAMBDA VALUES. THIS IS THE START OF THE N - R * C * ITERATION PROCESS. * C * * C **************************************************************** C DO 70,K=LOW$ATOM$SECT,HIGH$ATOM$SECT XYZP(K,4)=0.0D0 XYZP(K,5)=0.0D0 DO 71,J=1,3 71 XYZP(K,J)=AP1(K,J) 70 CONTINUE C C **************************************************************** C * * C * START CYCLE THROUGH THE MOLECULE. * C * * C **************************************************************** C DO 50,IC=1,NCYC C C **************************************************************** C * * C * START CYCLE THROUGH ATOMS. * C * * C **************************************************************** C DO 20,ITA=LOW$ATOM$SECT + 1,HIGH$ATOM$SECT-1 ITM1=ITA-1 ITP1=ITA+1 X1=XYZP(ITM1,1) X2=XYZP(ITA,1) X3=XYZP(ITP1,1) Y1=XYZP(ITM1,2) Y2=XYZP(ITA,2) Y3=XYZP(ITP1,2) Z1=XYZP(ITM1,3) Z2=XYZP(ITA,3) Z3=XYZP(ITP1,3) C C **************************************************************** C * * C * GET THE CURRENT BOND LENGTHS * C * * C **************************************************************** C BL1=DSQRT((X2-X1)**2 + (Y2-Y1)**2 + (Z2-Z1)**2) BL2=DSQRT((X3-X2)**2 + (Y3-Y2)**2 + (Z3-Z2)**2) C C **************************************************************** C * * C * TRAP FOR BAD DATA / BOND LENGTHS * C * * C **************************************************************** C IF ( BL2 .GE. 6.D0 ) THEN TYPE 600,BL1,ITA 600 FORMAT(1X,' YOU HAVE BAD DATA IN THE FORM OF A BOND LENGTH', 1 //' THE BOND LENGTH IS = ',D15.7,' THE ATOM SITE = ',I4) END IF C C **************************************************************** C * * C * GET THE FIRST DERIVATIVE (GRADIENT) AND JACOBIAN * C * * C **************************************************************** C CALL DER C C **************************************************************** C * * C * SET UP THE INVERSE AND DETERMINANT IN CASE ANOTHER * C * ALGORITHM FOR MATRIX INVERSION IS DESIRED. * C * * C **************************************************************** C CALL INV C C **************************************************************** C * * C * GET THE PARAMETER - CHANGE (SHIFT) VECTOR * C * * C **************************************************************** C CALL CHANGE C C **************************************************************** C * * C * ADD IN THE SHIFT TO GET THE NEW COORDINATES * C * * C **************************************************************** C XYZP(ITA,1)=XYZP(ITA,1) + C(1) XYZP(ITA,2)=XYZP(ITA,2) + C(2) XYZP(ITA,3)=XYZP(ITA,3) + C(3) XP=XYZP(ITA,1) YP=XYZP(ITA,2) ZP=XYZP(ITA,3) C C **************************************************************** C * * C * RESET THE LAMBDA VALUES FOR NEXT ITERATION * C * * C **************************************************************** C XYZP(ITA,4)= XYZP(ITA,4) +C(4) XYZP(ITA,5)= XYZP(ITA,5) +C(5) 20 CONTINUE C C **************************************************************** C * * C * GET THE BOND LENGTH AND POSITION DEVIATIONS * C * * C **************************************************************** C CALL DEV( &55,&50 ) 50 CONTINUE 55 CONTINUE C C **************************************************************** C * * C * OUTPUT THE INFORMATION FROM THIS SECTION OF THE * C * PROGRAM TO THE SPECIFIED FILE, I. E. AOP****.DAT * C * * C **************************************************************** C CALL SHIFT$OUT C C **************************************************************** C * * C * THIS IS THE BRANCH TO THE ANGLE SECTION OF THE * C * PROGRAM. * C * * C **************************************************************** C TYPE 700 700 FORMAT(1X,'DO YOU WANT TO PROCEED TO CALCULATE THE SHAPELY' 1 ' ANGLES (Y OR N)',$) ACCEPT 750,ANS 750 FORMAT(A1) IF( ANS .EQ. YES ) THEN TYPE 701,FILEOUTS 701 FORMAT(1X,' COORDINATES AND BOND LENGTHS IN FILE ',1A20 ) CALL ANGLES END IF 950 CONTINUE CLOSE(2) CLOSE(3) STOP END SUBROUTINE FILES C C **************************************************************** C * * C * SET UP THE VARIOUS FILES THAT ARE REQUIRED * C * BY BOTH SECTIONS OF THE PROGRAM. * C * * C **************************************************************** C C INCLUDE 'CHARCOM.FOR' C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C INCLUDE 'LBCOM1.FOR' C C **************************************************************** C * * C * GET THE INPUT FILE NAME FOR THE ALPHA - CARBON * C * ATOMS. THIS FILE IS CONSTRUCTED FROM THE * C * BROOKHAVEN DATA BANK FILES. * C * * C **************************************************************** C TYPE 200 200 FORMAT(1X,' ENTER PROTEIN ATOM INPUT,AP****.*** >',$) ACCEPT 205,FILEIN 205 FORMAT(A) C C **************************************************************** C * * C * GET THE FILE NAME FOR THE OUTPUT FROM THE SHIFT * C * SECTION OF THE PROGRAM. * C * * C **************************************************************** C TYPE 201 201 FORMAT(1X,' ENTER PROTEIN ATOM OUTPUT,AOP****.*** >',$) ACCEPT 205,FILEOUTS C C==>OPEN FILES C OPEN(UNIT=1,FILE=FILEIN,TYPE='OLD',READONLY) OPEN(UNIT=2,FILE=FILEOUTS,TYPE='NEW',CARRIAGECONTROL='LIST') TYPE 2003 2003 FORMAT(1X,' TITLE ON OUTPUT FILE (MAX. 80 CHAR.)',/' >',$) ACCEPT 205,TITLE C C **************************************************************** C * * C * GET THE OUTPUT FILE NAME FOR THE INFORMATION * C * FROM THE ANGLE SECTION OF THE PROGRAM. THIS * C * FILE WILL CONTAIN THE DATA NECESSARY TO CONSTRUCT * C * THE SHAPELY MODEL. * C * * C **************************************************************** C TYPE 1000 1000 FORMAT(1X,' ENTER OUTPUT FILE,BOP****.DAT >>',$) ACCEPT 205,FILEOUT C C==>OPEN THE FILES C OPEN(UNIT=4,FILE ='RESCOL.DAT',TYPE='OLD',READONLY) OPEN(UNIT=3,FILE =FILEOUT,TYPE='NEW',CARRIAGECONTROL='LIST') RETURN END SUBROUTINE ARRAY$SETUP C C **************************************************************** C * * C * SET UP THE INITIAL DATA STRUCTURES * C * * C **************************************************************** C C INCLUDE 'CHARCOM.FOR' C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C INCLUDE 'LBCOM1.FOR' C C C **************************************************************** C * * C * READ THE DATA. DATA MUST BE ALPHA CARBON * C * * C **************************************************************** C C C NA = 0 10 NA = NA + 1 READ(1, 211, END= 900) RES(NA), ORG$NUM(NA), (AP1(NA,J),J=1,3) GOTO 10 900 CONTINUE NA = NA - 1 211 FORMAT(1X,4X,A3,2X,I4,5X,3(D7.3,1X)) CLOSE(1) C C==>GET THE SEQUENCE AND COLOR CODES C TYPE 105 105 FORMAT(/,' THE CURRENT COLOR ASSIGNMENTS ARE:'/) DO 20,J=1,20 IRES$DSTRB(J) = 0 READ(4,103) RES1(J),RES3(J),COLOR(J) 20 TYPE 103,RES1(J),RES3(J),COLOR(J) CLOSE (4) IRES$TOTAL = 0 103 FORMAT(1X,A1,5X,A3,5X,A20,5X,I3,5X,I2) RETURN END SUBROUTINE SHIFT$OUT C INCLUDE 'CHARCOM.FOR' C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C C C **************************************************************** C * * C * THIS IS THE OUTPUT SECTION FROM THE SHIFT * C * SECTION OF THE PROGRAM * C * * C * * C **************************************************************** C WRITE(2,658) 658 FORMAT(1X,' OUTPUT FROM SHAPELY MODELS SHIFT PROGRAM ( LSM ) '/) WRITE(2,657) 657 FORMAT(1X,' (C) COPYRIGHT SHAPELY MODELS, 1981 '/) WRITE(2,650) ICH,AL,SSL,SSP,NUM$ATOM$SECT 650 FORMAT(1X,'CYCLE=',I3,2X,'AVERAGE LENGTH=',D15.7,2X, 1 ' SUM OF SQUARES OF LENGTH DEVIATIONS=',D15.7,/ 1 ' SUM OF SQUARES OF POSITION DEVIATIONS=',D15.7,2X, 1 ' NUMBER OF ATOMS=',I4/) WRITE(2,651)MAX$SHIFT,COR$FACT 651 FORMAT(1X,'MAX$SHIFT=',D15.7,5X,'COR$FACT=',D15.7/) WRITE(2,654) 654 FORMAT(1X,T4,'AA#2',T13,'NEW XYZ',T28,'NEW BOND LENGTH', 1 T55,'OLD XYZ',T69,'ORIGINAL BOND LENGTH',T96,'AA#1'/) DO 85,I=LOW$ATOM$SECT,HIGH$ATOM$SECT WRITE(2,655)I,(HOLD(I,J),J=1,3),NBLH(I),RES(I),(AP1(I,J),J=1,3), 1 ORBL(I),ORG$NUM(I) 655 FORMAT(1X,I4,2X,3(F5.1,1X),5X,F5.2,5X,A3,5X,3(F5.1,1X),2X,F5.2,20X, 1 I4) C C C **************************************************************** C * * C * * C * TEST FOR SHORT OR LONG NEW BOND LENGTHS * C * IF LAST POSITION,DO NOT CHECK * C * * C **************************************************************** C C IF (I .EQ. HIGH$ATOM$SECT ) GO TO 85 IF (NBLH(I) .LT. 3.7D0 .OR. NBLH(I) .GT. 3.9D0) THEN TYPE 656,I,NBLH(I) WRITE(2,656) I,NBLH(I) END IF 656 FORMAT(1X,'BAD BOND LENGTH...RES=',I3,2X,D15.7) 85 CONTINUE WRITE(2,700) 700 FORMAT(///) RETURN END SUBROUTINE DEV( *,* ) C C **************************************************************** C * * C * THIS ROUTINE GETS THE BOND LENGTH AND POSITION * C * DEVIATIONS AND ALSO TESTS FOR CONVERGENCE. * C * * C **************************************************************** C C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C C C **************************************************************** C * * C * * C * GET THE BOND LENGTH DEVIATIONS * C * * C **************************************************************** C BLEN=0.D0 BLEN1=0.D0 SSLDEV=0.D0 DO 80,IL=LOW$ATOM$SECT,HIGH$ATOM$SECT-1 BLEN1= DSQRT((XYZP(IL,1) - XYZP(IL+1,1))**2 1 + (XYZP(IL,2) - XYZP(IL+1,2))**2 1 + (XYZP(IL,3) - XYZP(IL+1,3))**2) NBL(IL)=BLEN1 SSLDEV=SSLDEV + (BLEN1 - TBL)**2 BLEN=BLEN + BLEN1 80 CONTINUE AVGLEN=BLEN / DFLOTJ(NUM$ATOM$SECT-1) SSLDEV=SSLDEV / DFLOTJ(NUM$ATOM$SECT-1) C C C **************************************************************** C * * C * * C * GET THE POSITION DEVIATIONS * C * * C **************************************************************** C POS=0.D0 DO 81,IP=LOW$ATOM$SECT + 1,HIGH$ATOM$SECT-1 81 POS= POS + DSQRT((AP1(IP,1) - XYZP(IP,1))**2 1 +(AP1(IP,2) - XYZP(IP,2))**2 1 +(AP1(IP,3) - XYZP(IP,3))**2) SSPOS=POS / DFLOTJ(NUM$ATOM$SECT-2) TYPE 650,IC,AVGLEN,SSLDEV,SSPOS,NUM$ATOM$SECT 650 FORMAT(1X,'CYCLE=',I3,2X,'AVERAGE LENGTH=',D15.7,2X, 1 ' SUM OF SQUARES OF LENGTH DEVIATIONS=',D15.7,/ 1 ' SUM OF SQUARES OF POSITION DEVIATIONS=',D15.7,2X, 1 ' NUMBER OF ATOMS=',I4/) C C C C **************************************************************** C * * C * * C * TEST FOR CONVERGENCE * C * * C **************************************************************** C IF (SSLDEV .GE. CONV) RETURN 1 CONV=SSLDEV AL=AVGLEN SSL=SSLDEV SSP=SSPOS ICH=IC DO 57,I=LOW$ATOM$SECT,HIGH$ATOM$SECT NBLH(I)=NBL(I) DO 57,J=1,3 57 HOLD(I,J)=XYZP(I,J) RETURN 2 END SUBROUTINE DER C C C **************************************************************** C * * C * * C * CALCULATES THE FIRST PARTIAL DERIVATIVES AND * C * THE ASSOCIATED JACOBIAN. * C * * C **************************************************************** C C C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL TWO=2.D0 LAM1=XYZP(ITA,4) LAM2=XYZP(ITA,5) C C C **************************************************************** C * * C * * C * CLEAR THE DERIVATIVE VECTOR AND THE JACOBIAN * C * MATRIX. * C * * C **************************************************************** C C DO 10,I=1,5 DF1(I)=0.D0 DO 10,J=1,5 10 JAC(I,J)=0.D0 C C C **************************************************************** C * * C * * C * FIRST DERIVATIVE SECTION: * C * DF1 IS THE GRADIENT VECTOR * C * * C **************************************************************** C C DF1(1)=( TWO * (X2-AP1(ITA,1)) + LAM1 * (X2-X1) / BL1 1 + LAM2 * (X3-X2) / BL2) DF1(2)=( TWO * (Y2-AP1(ITA,2)) + LAM1 * (Y2-Y1) / BL1 1 + LAM2 * (Y3-Y2) / BL2) DF1(3)=( TWO * (Z2-AP1(ITA,3)) + LAM1 * (Z2-Z1) / BL1 1 + LAM2 * (Z3-Z2) / BL2) DF1(4)=BL1-TBL DF1(5)=BL2-TBL C C C **************************************************************** C * * C * * C * JACOBIAN SECTION * C * * C **************************************************************** C C TL12= TWO * TBL + LAM1 + LAM2 JAC(1,1)=TL12 JAC(2,1)=0.0D0 JAC(3,1)=0.0D0 JAC(4,1)=(X2-X1) / BL1 JAC(5,1)=(X3-X2) / BL2 JAC(1,2)=0.0D0 JAC(2,2)=TL12 JAC(3,2)=0.0D0 JAC(4,2)=(Y2-Y1) / BL1 JAC(5,2)=(Y3-Y2) / BL2 JAC(1,3)=0.0D0 JAC(2,3)=0.0D0 JAC(3,3)=TL12 JAC(4,3)=(Z2-Z1) / BL1 JAC(5,3)=(Z3-Z2) / BL2 JAC(1,4)=(X2-X1) JAC(2,4)=(Y2-Y1) JAC(3,4)=(Z2-Z1) JAC(4,4)=0.0D0 JAC(5,4)=0.0D0 JAC(1,5)=(X3-X2) JAC(2,5)=(Y3-Y2) JAC(3,5)=(Z3-Z2) JAC(4,5)=0.0D0 JAC(5,5)=0.0D0 RETURN END SUBROUTINE INV C C C **************************************************************** C * * C * * C * SET UP THE INVERSE OF THE JACOBIAN MATRIX * C * IN CASE ANOTHER FORM OF MATRIX INVERSE ROUTINE * C * IS DESIRED. * C * * C **************************************************************** C C C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C C C **************************************************************** C * * C * * C * PLACE THE JACOBIAN MATRIX IN A HOLDING MATRIX * C * * C **************************************************************** C C DO 55,I=1,5 DO 55,J=1,5 55 AINV(I,J)=JAC(I,J) CALL MATINV RETURN END SUBROUTINE CHANGE C C C **************************************************************** C * * C * * C * THIS ROUTINE DETERMINES THE PARAMETER-CHANGE * C * VECTOR USING A JACOBIAN MATRIX. * C * * C **************************************************************** C C C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C C C **************************************************************** C * * C * * C * GET THE CHANGE VECTOR * C * * C **************************************************************** C C DO 10,I=1,5 C(I)=0.D0 DO 13,J=1,5 13 C(I)= (C(I) + (AINV(I,J) * ( - DF1(J)))) 10 CONTINUE C C C **************************************************************** C * * C * * C * CHECK FOR EXCESSIVE SHIFTS * C * NOTE: THIS IS WHERE THE TWO PARAMETERS MAY BE * C * TWEEKED,I.E. MAX$SHIFT AND COR$FACT * C * * C **************************************************************** C C DO 15,I=1,3 15 IF( DABS(C(I)) .GE. MAX$SHIFT ) C(I) = C(I) * COR$FACT RETURN END SUBROUTINE MATINV C C C **************************************************************** C * * C * * C * MATRIX INVERSION ROUTINE - GAUSS JORDAN METHOD * C * * C * ARGUMENTS: * C * AINV: INPUT MATRIX REPLACED BY ITS INVERSE * C * NORDER: ORDER OF THE MATRIX * C * DET: VALUE OF THE DETERMINANT * C * * C **************************************************************** C C REAL*8 AMAX,SAVE C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL DIMENSION IK(5),JK(5) NORDER=5 10 DET = 1.0D0 11 DO 100 K = 1,NORDER C C C **************************************************************** C * * C * * C * LOCATE LARGEST ELEMENT IN AINV * C * * C **************************************************************** C C AMAX = 0.0D0 21 DO 30 I=K,NORDER DO 30 J=K,NORDER 23 IF (DABS(AMAX) - DABS(AINV(I,J))) 24,24,30 24 AMAX = AINV(I,J) IK(K) = I JK(K) = J 30 CONTINUE C C C **************************************************************** C * * C * * C * INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN * C * AINV(K,K) * C * * C **************************************************************** C C 31 IF(AMAX) 41,32,41 32 DET = 0.0D0 GO TO 140 41 I = IK(K) IF (I-K) 21,51,43 43 DO 50 J=1,NORDER SAVE = AINV(K,J) AINV(K,J) = AINV(I,J) 50 AINV(I,J) = -SAVE 51 J = JK(K) IF (J-K) 21,61,53 53 DO 60 I=1,NORDER SAVE = AINV(I,K) AINV(I,K) = AINV(I,J) 60 AINV(I,J) = -SAVE C C C **************************************************************** C * * C * * C * ACCUMULATE ELEMENTS OF INVERSE MATRIX * C * * C **************************************************************** C C 61 DO 70 I=1,NORDER IF(I-K) 63,70,63 63 AINV(I,K) = -AINV(I,K)/AMAX 70 CONTINUE 71 DO 80 I=1,NORDER DO 80 J=1,NORDER IF(I-K) 74,80,74 74 IF(J-K) 75,80,75 75 AINV(I,J) = AINV(I,J) + AINV(I,K)*AINV(K,J) 80 CONTINUE 81 DO 90 J=1,NORDER IF(J-K) 83,90,83 83 AINV(K,J) = AINV(K,J)/AMAX 90 CONTINUE AINV(K,K) = 1.0/AMAX 100 DET = DET*AMAX C C C **************************************************************** C * * C * * C * RESTORE ORDERING OF MATRIX * C * * C **************************************************************** C C 101 DO 130 L=1,NORDER K = NORDER - L + 1 J = IK(K) IF(J-K) 111,111,105 105 DO 110 I=1,NORDER SAVE = AINV(I,K) AINV(I,K) = -AINV(I,J) 110 AINV(I,J) = SAVE 111 I = JK(K) IF(I-K) 130,130,113 113 DO 120 J=1,NORDER SAVE = AINV(K,J) AINV(K,J) = -AINV(I,J) 120 AINV(I,J) = SAVE 130 CONTINUE 140 CONTINUE RETURN END SUBROUTINE ANGLES C C **************************************************************** C * * C * ROUTINE TO CALCULATE THE TWO ANGLES ( BEND AND * C * TORSIONAL ) NEEDED TO CONSTRUCT ALPHA - CARBON * C * MODELS USING SHAPELY MODEL COMPONENTS. THE * C * DATA COME FROM THE SHIFT SECTION OF THE PROGRAM. * C * * C * NOTE: THE BEND ANGLE IS LISTED AS THETA AND THE * C * TORSIONAL ANGLE AS TAU. * C **************************************************************** C C C INCLUDE 'CHARCOM.FOR' C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C INCLUDE 'LBCOM1.FOR' IRES$SECT$TOTAL = 0 DO 20,I=1,20 20 IRES$DSTRB(I) = 0 C C C **************************************************************** C * * C * * C * GET THE ALPHA CARBON COORDINATES AND RESIDUE TYPE * C * * C **************************************************************** C C C DO 11,I=LOW$ATOM$SECT,HIGH$ATOM$SECT PRES3(I) = RES(I) DO 11,J=1,3 XYZ(I,J) = XYZP(I,J) 11 CONTINUE C C C **************************************************************** C * * C * * C * MATCH COLOR CODES,I.E. PROTEIN RESIDU VS. SHAPELY * C * ASSIGNMENTS. * C * * C **************************************************************** C C DO 12,I=LOW$ATOM$SECT,HIGH$ATOM$SECT DO 13,J=1,20 IF( PRES3(I) .NE. RES3(J)) GOTO 13 IRES$SECT$TOTAL = IRES$SECT$TOTAL + 1 IRES$DSTRB$TOTAL(J) = IRES$DSTRB$TOTAL(J) + 1 IRES$DSTRB(J)=IRES$DSTRB(J) + 1 PRES1(I) = RES1(J) PCOL(I) = COLOR(J) IND$COL(I) = J GOTO 12 13 CONTINUE 12 CONTINUE C C C **************************************************************** C * * C * * C * TRAP FOR LOW ATOM COUNT IN SECTION * C * ALSO SET THE VALUES FOR START OF LOOP * C * * C **************************************************************** C C IF( LOW$ATOM$FLAG .EQ. 1 ) RETURN C LA1 = LOW$ATOM$SECT LA2 = LOW$ATOM$SECT + 1 LA3 = LOW$ATOM$SECT + 2 DO 15,J=1,3 VPREV(J)=XYZ(LA2,J)-XYZ(LA1,J) 15 VNOW(J)=XYZ(LA3,J)-XYZ(LA2,J) C C==>DISTANCE OF PREVIOUS VECTOR C DSTPRV=SQRT(VPREV(1)**2 + VPREV(2)**2 + VPREV(3)**2) C C==>GET THE CROSS PRODUCT C CALL CROSS(VPREV,VNOW,CRSPRV) C C==>DISTANCE OF THE CROSS PRODUCT C DSTCRP=SQRT(CRSPRV(1)**2 + CRSPRV(2)**2 + CRSPRV(3)**2) C C==>START THE CYCLE THROUGH THE MOLECULE C DO 50,N=LOW$ATOM$SECT + 1,HIGH$ATOM$SECT-1 C C==>CALCULATE THE BEND ANGLE C DSTNOW=SQRT(VNOW(1)**2 + VNOW(2)**2 + VNOW(3)**2) C C==>GET THE VECTOR DOT VECTOR PRODUCT C VDOTV=VPREV(1)*VNOW(1) + VPREV(2)*VNOW(2) + VPREV(3)*VNOW(3) C C==>GET THE BEND ANGLE C IBEND=ABS(57.2958 * (ACOS(ABS(VDOTV/(DSTPRV*DSTNOW))))) + 0.5 IF( VDOTV .LT. 0.0) IBEND = 180 - IBEND C C==>CALCULATE THE DIHEDRAL ANGLE,CHECK HANDEDNESS C DO 25,J=1,3 25 VNEXT(J)=XYZ(N+2,J) - XYZ(N+1,J) CALL CROSS(VNOW,VNEXT,CRSNXT) DSTCRN=SQRT(CRSNXT(1)**2 + CRSNXT(2)**2 + CRSNXT(3)**2) CDOTC=CRSPRV(1)*CRSNXT(1) + CRSPRV(2)*CRSNXT(2) 1 + CRSPRV(3)*CRSNXT(3) C C==>ON OCCASION THE ARGUMENT FOR ACOS IS SLIGTLY GREATER THAN ONE C==>E.G. 1.00000002, WHICH CAUSES AN ERROR TRAP FOR THIS BELOW: C ARG= ABS( CDOTC / (DSTCRP * DSTCRN) ) IF ( ARG .GT. 1.001 ) STOP 'BAD ANGLE IN ANGLES' IF ( ARG .GT. 1.0000000 ) ARG= 1.0000000 DIHEDR = ABS( 57.2958 * ACOS( ARG )) IF( CDOTC .LT.0.0) DIHEDR = 180. - DIHEDR CALL CROSS(CRSPRV,CRSNXT,CCRSC) VDOTCC=VNOW(1)*CCRSC(1) + VNOW(2)*CCRSC(2) 1 + VNOW(3)*CCRSC(3) IF( VDOTCC .LT. 0.0) DIHEDR = 360. - DIHEDR IF( DIHEDR .GE. 360.0) DIHEDR = DIHEDR - 360.0 C C==>SAVE ANGLES C IDIH(N)=DIHEDR+ 0.5 IBEN(N)=IBEND C C==>SET THE VALUES FOR THE NEXT PASS C DIHPRV=DIHEDR DO 55,J=1,3 VPREV(J)=VNOW(J) VNOW(J)=VNEXT(J) 55 CRSPRV(J)=CRSNXT(J) DSTPRV=DSTNOW DSTCRP=DSTCRN 50 CONTINUE CALL ANGLE$OUT RETURN END SUBROUTINE ANGLE$OUT C C C **************************************************************** C * * C * * C * PRODUCES THE OUTPUT FROM THE ANGLE SECTION * C * OF THE PROGRAM. THIS INFORMATION IS USED * C * IN THE CONSTRUCTION OF THE SHAPELY MODEL. * C * * C **************************************************************** C C C INCLUDE 'CHARCOM.FOR' C INCLUDE 'LSCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL C INCLUDE 'LBCOM1.FOR' C C==>OUTPUT SECTION C WRITE (3,2010) 2010 FORMAT (///, ' ******************** SHAPELY', # ' MODELS ANGLE CALCULATION *****************',//, 1 ' (C) COPYRIGHT SHAPELY MODELS,1982 ',//, # ' NOTE THAT THE N-TERMINAL SIDE IS THE MALE CONNECTOR. THE', # ' BEND ANGLE IS AT THE ALPHA CARBON ATOM. '/ # ' THE TORSIONAL ANGLE IS FOLLOWING THE ALPHA CARBON.' # /, ' THE ANGLE LIST IS BEND ANGLE ( THETA ) FOLLOWED BY ', # 'TORSIONAL ANGLE ( TAU ) ') WRITE (3,2005) TITLE 2005 FORMAT (//, ' TITLE- ', A80,// ) WRITE(3,107) 107 FORMAT(/,T1,'RES NO.',T9,'TYPE(3)',T17,'TYPE(1)',T25,'THETA', 1 T35,'TAU',T47,'COLOR',T66,'COLOR NO.',T91,'AA#1',/) IBREAK = 0 DO 60,I=LOW$ATOM$SECT,HIGH$ATOM$SECT IBREAK = IBREAK + 1 IF( IBREAK .EQ. 10 ) THEN WRITE(3,301) I,PRES3(I),PRES1(I),IBEN(I),IDIH(I),PCOL(I),IND$COL(I), 1 ORG$NUM(I) 301 FORMAT(1X,I4,'*',4X,A3,5X,A1,2(5X,I4),5X,A20,5X,I2,20X,I4) IBREAK = 0 GOTO 60 END IF WRITE(3,300) I,PRES3(I),PRES1(I),IBEN(I),IDIH(I),PCOL(I), 1 IND$COL(I),ORG$NUM(I) CALL DIST (I ) 60 CONTINUE 300 FORMAT(1X,I4,5X,A3,5X,A1,2(5X,I4),5X,A20,5X,I2,20X,I4) WRITE(3,104) 104 FORMAT(/,' SHAPELY MODELS CURRENT(RESCOL.DAT) COLOR ASSIGNMENTS'/, 1 ' AND AMINO ACID/COLOR DISTRIBUTION') DO 26,J=1,20 IRES$TOTAL = IRES$TOTAL + IRES$DSTRB(J) 26 WRITE(3,103)RES1(J),RES3(J),COLOR(J),IRES$DSTRB(J),J WRITE(3,106) IRES$SECT$TOTAL 106 FORMAT(/' THE TOTAL NUMBER OF RESIDUES IN THE SECTION =',I4) 103 FORMAT(1X,A1,5X,A3,5X,A20,5X,I3,5X,I2) WRITE(3,700) 700 FORMAT(///) IF ( IJK .EQ. SECT$NUM ) THEN WRITE(3,700) DO 27, I=1,20 WRITE(3,103) RES1(I),RES3(I),COLOR(I),IRES$DSTRB$TOTAL(I),I 27 CONTINUE WRITE(3,108) IRES$TOTAL 108 FORMAT(/' THE TOTAL NUMBER OF RESIDUES IN THE MOLECULE = ',I4) END IF RETURN END SUBROUTINE CROSS(A,B,C) DIMENSION A(3),B(3),C(3) C(1) = A(2) * B(3) - A(3) * B(2) C(2) = A(3) * B(1) - A(1) * B(3) C(3) = A(1) * B(2) - A(2) * B(1) RETURN END SUBROUTINE DIST(IM) C C C **************************************************************** C * * C * * C * CALCULATES THE DISTANCE TO NEIGHBORS OF ATOM I, * C * THAT ARE WITHIN A SPECIFIED RADIUS. * C * * C **************************************************************** C C C INCLUDE 'CHARCOM.FOR' C INCLUDE 'LBCOM1.FOR' CHARACTER*1 YES,ANS,RES1,PRES1,IW CHARACTER*3 RES3,PRES3,RES CHARACTER*20 COLOR,PCOL CHARACTER*80 TITLE CHARACTER*10 FILEIN CHARACTER*11 FILEOUT,FILEOUTS REAL*8 DF1,JAC,AINV,C,DET,TBL,AP1,HOLD,NBLH,XYZP,ORBL, 1 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,TL12,MAX$SHIFT,COR$FACT, 1 BLEN,POS,AVGLEN,SSPOS,SSLDEV,BLEN1,NBL,CONV, 1 AL,SSL,SSP,LAM1,LAM2 INTEGER*4 ORG$NUM,SECT,SECT$NUM,HIGH$ATOM$SECT COMMON/BLK1/ NA,DF1(5),JAC(5,5),AINV(5,5),C(5),DET,TBL,LAM1,LAM2 COMMON/BLK2/AP1(1000,3),XYZP(1000,5),ORBL(1000) COMMON/BLK3/X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,BL1,BL2,ITA,MAX$SHIFT,COR$FACT COMMON/BLK8/ORG$NUM(1000),SECT(10,4),NUM$ATOM$SECT,SECT$NUM COMMON/BLK9/RES(1000) COMMON/BLK10/NBL(1000),HOLD(1000,3),NBLH(1000) COMMON/BLK12/CONV,AL,SSL,SSP COMMON/BLK13/ICH,IC,IJK,LOW$ATOM$FLAG COMMON/BLK14/AVGLEN,SSLDEV,SSPOS COMMON/BLK16/FILEOUTS COMMON/BLK17/TITLE COMMON/BLK15/HIGH$ATOM$SECT,LOW$ATOM$SECT COMMON/BLK4/RES1(20),RES3(20),COLOR(20) COMMON/BLK5/PRES1(1000),PRES3(1000),PCOL(1000) COMMON/BLK6/IBEN(1000),IDIH(1000),XYZ(1000,3) COMMON/BLK7/VPREV(3),VNOW(3),VNEXT(3),CCRSC(3),CRSNXT(3),CRSPRV(3) COMMON/BLK11/IRES$DSTRB(20),IND$COL(1000),IRES$DSTRB$TOTAL(20) COMMON/BLK18/IRES$TOTAL,IRES$SECT$TOTAL DIMENSION IH(3),IW(3) DATA IW/' ','>',' '/ C C C **************************************************************** C * * C * * C * NOTE: THIS VERSION LOOKS UP TO 10 AMINO ACID * C * UNITS BACK. THAT IS THE CHAIN IS SCANNED * C * FROM THE N - TERMINUS TO THE I-10 TH. UNIT * C * * C **************************************************************** C C IM5 = IM - 10 DO 300,J=1,IM5 DS = 0 DO 301,K=1,3 D = XYZ(IM,K) - XYZ(J,K) D = D * D DS = D + DS 301 CONTINUE DS = SQRT( DS ) IF( DS .GT. 5.0 .AND. DS .LT. 7.0 ) GOTO 304 GOTO 300 304 IH(1) = J - 1 IH(3) = J + 1 IH(2) = J IF( IH(1) .EQ. 0 ) GOTO 300 WRITE(3,3000) DS 3000 FORMAT(/T70,' THE NEIGHBOR AT A',F5.1,' CM. RADIUS IS:') DO 50,I=1,3 WRITE(3,3001) IW(I),IH(I),PRES3(IH(I)),PRES1(IH(I)),PCOL(IH(I)) 50 CONTINUE 3001 FORMAT(T70,A1,I3,5X,A3,3X,A1,2X,A20) 300 CONTINUE RETURN END