C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C C ============== PROGRAM VOLUMEX C ============== C C (CORRECTED TO 14 APR 83) C C (The option to calculate cavity volumes is not activated in this C version of the program.) C C To calculate the volume of a series of specified atoms in a list of C coordinates by the voronoi, or various modified voronoi, procedures. C C The input file for this program is assumed to have been prepared C through the program access with or without actual calculation of C the accessible areas of the atoms. the key array in this file C indicates the atoms to be excluded or included in the calculation C and which included atoms are to have their volumes calculated. C C Functions or subroutines called: getnum,block C shell, filbox, atmpol, polvol C (with d_lines also prntbx,prntpn) C (note that, if wanted, the image of icl is C prepared in shell with d_lines) C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM,ITEMP CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) COMMON/FMRTST/NCALLS(25),JLEVEL,IMSG,IERROR(2,5),ICHK,JCHK CHARACTER*6 KEYWRD EXTERNAL VOLBLK LOGICAL VERBOSE,LOGSHELL,LOGWATER INTEGER NPARM PARAMETER (NPARM=200) COMMON /LOGCOM/ VERBOSE REAL FVALUE(NPARM) INTEGER IBEG(NPARM),IEND(NPARM), + ITYP(NPARM),IDEC(NPARM) LOGICAL LEND,LOGCALC CHARACTER CVALUE(NPARM)*4,LINE*255,KEYP*4 C C---- default values C CE = 2.80 ARANGE = 6.50 RDELTA = 0.0 C C---- assign input and output file names C logical unit #1 - input file C logical unit #2 - output file : notes and data C logical unit #3 - output file : error messages C logical unit #4 - input/output file : shell coords. C VERBOSE = .FALSE. LOGSHELL = .FALSE. LOGCALC = .TRUE. LOGWATER = .FALSE. METHOD = 1 LDUMX = 80 IFAILX = 0 C C ************************************************* CALL CCPFYP CALL CCPRCS(6,'VOLUME','$Date: 2002/08/12 09:01:38 $') C ************************************************** C 6655 CONTINUE NTOK=NPARM LINE = ' ' C C *********************************************************** CALL PARSER(KEYP,LINE,IBEG,IEND,ITYP, +FVALUE,CVALUE,IDEC,NTOK,LEND, + .TRUE.) C *********************************************************** C IF (LEND .OR. KEYP.EQ.'END') GOTO 8877 C IF (KEYP.EQ.'EDGE') THEN IF (NTOK.NE.2 .OR. ITYP(2).NE.2) + CALL CCPERR (1, 'EDGE needs single numeric argument') CE = FVALUE(2) ELSE IF (KEYP.EQ.'SRAD') THEN IF (NTOK.NE.2 .OR. ITYP(2).NE.2) + CALL CCPERR (1, 'SRAD needs single numeric argument') ARANGE = FVALUE(2) ELSE IF (KEYP.EQ.'METH') THEN IF (NTOK.NE.2) CALL CCPERR (1, 'Missing method') KEYP = CVALUE(2) CALL CCPUPC(KEYP) IF (KEYP.EQ.'VORO') THEN METHOD = 1 ELSE IF (KEYP.EQ.'RICH') THEN METHOD = 2 ELSE IF (KEYP.EQ.'RADI') THEN METHOD = 3 ELSE CALL CCPERR (1, 'Invalid method: '//KEYP) END IF ELSE IF (KEYP.EQ.'OUTE') THEN IF (NTOK.NE.2 .OR. ITYP(2).NE.2) + CALL CCPERR (1, 'OUTER needs single numeric argument') RDELTA = FVALUE(2) ELSE IF (KEYP.EQ.'SHEL') THEN LOGSHELL = .TRUE. ELSE IF (KEYP.EQ.'NOSH') THEN LOGSHELL = .FALSE. ELSE IF (KEYP.EQ.'NOCA') THEN LOGCALC = .FALSE. ELSE IF (KEYP.EQ.'CALC') THEN LOGCALC = .TRUE. ELSE IF (KEYP.EQ.'VERB') THEN VERBOSE = .TRUE. ELSE IF (KEYP.EQ.'NOVE') THEN VERBOSE = .FALSE. ELSE IF (KEYP.EQ.'WATE') THEN LOGWATER = .TRUE. ELSE IF (KEYP.EQ.'NOWA') THEN LOGWATER = .FALSE. ELSE IF (KEYP.EQ.'RUN') THEN GO TO 8877 END IF C GOTO 6655 8877 CONTINUE WRITE(6,*) ' METHOD CHOSEN IS ',METHOD C C---- pick search radius,arange. C C---- select shell parameter, ce C WRITE(6,*) + ' Current value for length of the edge of the "shell"' WRITE(6,5002) CE 5002 FORMAT(' Lattice, CE, is : ',F5.2) C C---- pick search radius,arange. C WRITE(6, 5010) ARANGE 5010 FORMAT(' Current value of Search Radius, ARANGE, is: ',F5.2) C C---- decide method to be used for volume calculation C WRITE(6,*)' 1) Standard Voronoi Procedure' WRITE(6,*)' 2) Richards" Method B' WRITE(6,*)' 3) Radical Plane Procedure' C C---- pick value to position outer surface, rdelta C WRITE(6, 5020) RDELTA 5020 FORMAT(' For Surface Volume Adjustment,',/, 1 ' Current Value of RDELTA is: ',F5.2) C C---- Open files C C ************************************************ CALL CCPDPN(1,'XYZIN','READONLY','F',LDUMX,IFAILX) CALL CCPDPN(2,'XYZOUT','NEW','F',LDUMX,IFAILX) CALL CCPDPN(3,'ERROUT','NEW','F',LDUMX,IFAILX) CALL CCPDPN(4,'SHELLFILE','UNKNOWN','F',LDUMX,IFAILX) C ************************************************ C WRITE(2,5106) METHOD 5106 FORMAT(/,' The METHOD of Calculation is: ',I2) C C---- initialize error variables C JLEVEL=1 C DO 90 J=JLEVEL,5 DO 90 I=1,2 90 IERROR(I,J)=0 C ICHK=0 IMSG=1 IERROR(1,JLEVEL)=19 C C---- load common block /coords/ from input file C 100 READ(1,5001,END=5200) KEYWRD 5001 FORMAT(A6,74A1) CALL CCPUPC(KEYWRD) C IF(KEYWRD.NE.' BEGIN') GO TO 100 GO TO 150 C C---- error write C 5200 CALL CCPERR(1,' error in file format. no "begin" found') C 150 LSTATM=0 LAST=0 NEXT=0 C DO 200 I=1,MAXATM READ(1,5000,END=300) KEY(I),NUMCHN(I),NUMFIL(I), 1 ATM(I),RES3(I),SEQNUM(I),X(I),Y(I),Z(I),RVDW(I), 2 RCOV(I) 5000 FORMAT(1X,I2,7X,I2,1X,I2,1X,2A4,1X,I3,3F8.3, 1 2F5.2,2F6.1,F5.2) C C---- Check for and reset seqnum to increase constantly from C beginning of file to end. bnl files have some odd entries C for heme,hoh,etc. C NOW=SEQNUM(I) C IF(NOW.EQ.LAST) THEN SEQNUM(I)=NEXT ELSE IF(NOW.GT.LAST.OR.NOW.LT.LAST) THEN NEXT=NEXT+1 SEQNUM(I)=NEXT LAST=NOW END IF C LSTATM=LSTATM+1 200 CONTINUE C C---- all arrays have now been loaded from input file C determine min and max values of coord list. C 300 LSTRES=SEQNUM(LSTATM) XMIN=9999 YMIN=9999 ZMIN=9999 XMAX=-9999 YMAX=-9999 ZMAX=-9999 C DO 350 I=1,LSTATM IF(X(I).LT.XMIN) XMIN=X(I) IF(Y(I).LT.YMIN) YMIN=Y(I) IF(Z(I).LT.ZMIN) ZMIN=Z(I) IF(X(I).GT.XMAX) XMAX=X(I) IF(Y(I).GT.YMAX) YMAX=Y(I) IF(Z(I).GT.ZMAX) ZMAX=Z(I) 350 CONTINUE C WRITE(2,5104) LSTATM 5104 FORMAT(/,' Total number of protein atoms in input list is: ',I5) NPAT = LSTATM WRITE(6,*)' Input data loading complete' C C---- the protein is considered to consist of all atoms flagged with C key(i) = 0 or 1. (i.e. all atoms in the list flaged with -1 C are simply deleted completely from shell and all subsequent C calculations in this run. C C---- Set up hypothetical 'solvent' shell to define the protein C surface. the cubic lattice used for this purpose has an C edge length of ce. C 400 CE2=CE/2. C IF (LOGSHELL) CALL SHELL C C---- read in shell coordinates (new or old) from the file on lun 4 C IF (.NOT. LOGWATER) THEN MAXSOL = 0 NSOL = 0 NMIN2 = 0 NTOT = 0 ELSE READ(4,9001) MAXSOL, NSOL, NMIN2, NTOT 9001 FORMAT(4I5) C DO 425 J=LSTATM+1, LSTATM+NSOL READ(4,5000) KEY(J),NUMCHN(J),NUMFIL(J),ATM(J),RES3(J), 1 SEQNUM(J),X(J),Y(J),Z(J),RVDW(J),RCOV(J) 425 CONTINUE END IF C C---- cavity ("-2") cubes not included in this read as given C WRITE(6,*) ' shell loading complete' C IF (.NOT.LOGCALC) THEN CLOSE(UNIT=1,STATUS='KEEP') CLOSE(UNIT=2,STATUS='KEEP') CLOSE(UNIT=3,STATUS='KEEP') CLOSE(UNIT=4,STATUS='KEEP') CALL CCPERR(0, 'Normal termination') ENDIF C C---- start volume calculation for all atoms flagged with key(i)=0 C WRITE(2,5105) 5105 FORMAT(/,' # SER# CH FI ATM RES RES#',4X,'X',7X,'Y',7X,'Z', 1 3X,' RVDW RCOV VOLUME NSOL',/) WRITE(2,5107) 5107 FORMAT(' BEGIN') NCALC=0 TOTVOL=0 C DO 1000 I=1,LSTATM IF(KEY(I).NE.0) GO TO 1000 C C---- load ao(3) with the coordinates of the atom whose volume is C to be calculated (i.e. the "origin" atom.). C also load other parameters of ao. C AO(1)=X(I) AO(2)=Y(I) AO(3)=Z(I) AOVDW=RVDW(I) AOCOV=RCOV(I) AOSEQ=SEQNUM(I) IF(METHOD.NE.2) GOTO 450 C be careful about (updated) AONUM as it's in common ITEMP = AONUM C ******************** CALL GETNUMX(I,NRES,ITEMP) C ******************** AONUM = ITEMP C C---- error write C IF(AONUM.EQ.0.AND.METHOD.EQ.2) THEN WRITE(3,5108) 5108 FORMAT(' VOLUME: AONUM=0') GO TO 620 END IF C C---- initialize the arrays box,ibox,polyhn,ipolyh,vervol C 450 IMSG=1 ICHK=0 NASPHR=0 C DO 889 JDO=1,100 IBOX(1,JDO) = 0 IBOX(2,JDO) = 0 VERVOL(JDO) = 0.0 BOX(1,JDO) = 0.0 BOX(2,JDO) = 0.0 BOX(3,JDO) = 0.0 BOX(4,JDO) = 0.0 BOX(5,JDO) = 0.0 IPOLYH(1,JDO) = 0 IPOLYH(2,JDO) = 0 IPOLYH(3,JDO) = 0 IPOLYH(4,JDO) = 0 POLYHN(1,JDO) = 0.0 POLYHN(2,JDO) = 0.0 POLYHN(3,JDO) = 0.0 889 CONTINUE C 500 CALL FILBOX IF(ICHK.NE.-1) GO TO 520 C C---- error write C WRITE(3,5005) I 5005 FORMAT(' Box capacity exceeded in filbox for atom no.',i10) GO TO 620 520 CONTINUE ID=I CALL PRNTBX(ID) C JVMAX=100 IERROR(2,JLEVEL)=600 C 600 CALL ATMPOL(JVMAX) C IF(IMSG.EQ.6) GO TO 900 IF(IMSG.EQ.7) GO TO 900 IF(IMSG.EQ.8) GO TO 900 IF(ICHK.NE.-1) GO TO 650 C C---- error write C 620 WRITE(3,5006) I 5006 FORMAT(' CALCULATION OMITTED FOR ATOM NO.',I10,' GO TO NEXT ', 1 'ATOM IN INPUT LIST.',/) GO TO 900 650 CONTINUE C CALL PRNTPN(JVMAX) C C---- calculate the volume of the limiting polyhedron C 700 CALL POLVOL(JVMAX) C NCALC=NCALC + 1 TOTVOL=TOTVOL+VOLUME C C---- determine the number of solvent molecules C used in forming the polyhedron around ao. C NS=0 DO 750 JJ=1,100 IF(.NOT.VRTATM(JJ)) GO TO 750 IF(IBOX(2,JJ).EQ.1) NS=NS+1 750 CONTINUE C C---- write to output file for each atom. C IF(MOD(NCALC,100).EQ.0) WRITE(6, 5008) NCALC 5008 FORMAT(1X,I5,' ATOMS COMPLETED') GO TO 950 C C---- (note that if program error occurs C earlier, volume and ns are both set C equal to 1000. in normal calculation actual C values are entered.) C 900 VOLUME = 1000. NS = 1000 950 WRITE (2,5007) NCALC,I,NUMCHN(I),NUMFIL(I),ATM(I),RES3(I), 1 SEQNUM(I),X(I),Y(I),Z(I),RVDW(I),RCOV(I),VOLUME,NS 5007 FORMAT (2I5,2I3,2(A4,2X),I5,3F8.3,2F5.2,F8.3,I5) C 1000 CONTINUE C WRITE(2,5030) TOTVOL 5030 FORMAT(' TOTAL VOLUME OF CALCULATED ATOMS =',F10.2) C CLOSE(UNIT=1,STATUS='KEEP') CLOSE(UNIT=2,STATUS='KEEP') CLOSE(UNIT=3,STATUS='KEEP') CLOSE(UNIT=4,STATUS='KEEP') CALL CCPERR(0, 'Normal termination') END C C C C ================= BLOCK DATA VOLBLK C ================= C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) C C---- the array indx1(i) provides a pointer to the location in C array atnama where the data for residue i start. C SAVE DATA INDX1/ 1, 6, 10, 15, 20, 28, 33, 39, 45, 53, 56, 60, 1 69, 74, 85, 94, 105, 107, 110, 123, 127, 132, 134, 136, 2 138,140,185, 3*0/ C C---- the array indx2(i) provides a pointer to the location in C array ibnd1(1) where the data for residue i start. C DATA INDX2/1,6,10,15,20,29,34,40,46,54,57,61,71,76,89,98, 1 109,111,114,129,134,140,190,7*0/ C C---- for the arrays ibnd1() and atnama(): C note that for asn,gln,his 2 side chain lists are given C first for the 'a'designation,second for the complete n,c,o C designation. note that increment number will only give the C 'usual' atom number for the first list in each case. C C---- in the array ibnd1(i) the first entry for each residue C gives the increment in the index i for the last entry C for that residue. other entries give bond codes for each C bond in the main or side chains. C bond code=residue difference(0 or 1)*10 000+ lower numbered C atom in bond pair * 100 + higher numbered atom in bond pair. C DATA IBND1/4,102,103,304,10203,3,105,506,507,4,105,506,607,608, 1 4,105,506,507,608,8,105,506,607,608,709,810,911,1011, 2 4,105,506,607,608,5,105,506,607,708,709,5,105,506,607,708,809, 3 7,105,506,607,708,809,910,911,2,105,506,3,105,506,507, 4 9,105,506,607,608,709,810,911,1011,1112,4,105,506,607,708, 5 12,105,506,607,608,709,810,811,910,1012,1113,1214,1314, 6 8,105,506,607,608,109,910,1011,1012, 7 10,105,506,607,708,709,110,1011,1112,1213,1214,1,105,2,105,506, 8 14,105,506,607,608,709,810,910,111,1112,1213,1214,1315,1416,1516, 9 4,105,207,506,607,5,105,207,506,607,608,49,106,117,125,133,207, 1 237,310,318,421,426,529,534,607,610,708,809,812,910,911,1213, 2 1314,1415,1416,1718,1721,1819,1920,1922,2021,2023,2324,2526,2529, 3 2627,2728,2730,2829,2831,3132,3334,3337,3435,3536,3538,3637,3940, 44041,4142,4143,2,103,304,8*0/ C C---- the array restyp(i) contains the three letter name of the C residue i. C DATA RESTYP/'GLY ','VAL ','LEU ','ILE ','PHE ','ASP ','GLU ', 1 'LYS ','ARG ','SER ','THR ','TYR ','MET ','TRP ','ASN ','GLN ', 2 'ALA ','CYS ','HIS ','PRO ','HYP ','EXT ','DUM ','WAT ', 3 'HOH ','HEM ','ACE ', 3*' '/ C C---- in the array atnama(i) the first entry for each residue C gives the increment in the index i for the last entry C for that residue. other entries give the atom names in C 1-4 alphameric charachters.(see description of ityp1) C DATA KATNAMA/4, 4*0,3,3*0,4,4*0,4,4*0,7,7*0,4,4*0, + 5, 5*0, 5, 5*0, 7,7*0, 2, 2*0, 3, 3*0, 8,8*0, + 4, 4*0, 10,10*0,8,8*0,10,10*0,1,0,2,0,0, + 12,12*0, 3, 3*0,4,4*0,1,0,1,0,1,0,1,0,43,43*0, + 1,0,15*0/ DATA (ATNAMA(I),I=1,109)/'4. ', +'CA ', 'N ', 'C ', 'O ', '3. ', 'CB ', + 'CG1 ', 'CG2 ', '4. ', 'CB ', 'CG ', 'CD1 ', 'CD2 ', + '4. ', + 'CB ', 'CG1 ', 'CG2 ', 'CD1 ', '7. ', 'CB ', 'CG ', + 'CD1 ', 'CD2 ', 'CE1 ', 'CE2 ', 'CZ ', '4. ', 'CB ', + 'CG ', 'OD1 ', 'OD2 ', '5. ', 'CB ', 'CG ', 'CD ', + 'OE1 ', 'OE2 ', '5. ', 'CB ', 'CG ','CD ', 'CE ', + 'NZ ', '7. ', 'CB ', 'CG ', 'CD ', 'NE ', 'CZ ', + 'NH1 ', 'NH2 ', '2. ', 'CB ', 'OG ', + '3. ', 'CB ', 'OG1 ', + 'CG2 ', '8. ', 'CB ', 'CG ', 'CD1 ', 'CD2 ', 'CE1 ', + 'CE2 ', 'CZ ', 'OH ', '4. ', 'CB ', 'CG ', 'SD ', + 'CE ', '10. ', 'CB ', 'CG ', 'CD1 ', 'CD2 ', 'NE1 ', + 'CE2 ', 'CE3 ', 'CZ2 ', 'CZ3 ', 'CH2 ', + '8. ', 'CB ', 'CG ', + 'AD1 ', 'AD2 ', 'CB ', 'CG ', 'OD1 ', 'ND2 ', '10. ', + 'CB ', 'CG ', 'CD ', 'AE1 ', 'AE2 ', 'CB ', 'CG ', + 'CD ', 'OE1 ', 'NE2 ', '1. ', 'CB ', ' 2.', + 'CB ', 'SG '/ DATA (ATNAMA(I),I=110,200)/ + '12. ', 'CB ', 'CG ', 'AD1 ', 'AD2 ', 'AE1 ', 'AE2 ', + 'CB ', 'CG ', 'ND1 ', 'CD2 ', 'CE1 ', 'NE2 ', '3. ', + 'CB ', 'CG ', 'CD ', '4. ', 'CB ', 'CG ', 'CD ', + 'OD ', '1. ', 'EXT ', '1. ', 'DUM ', '1. ', +'WAT ', '1. ', 'HOH ', + '43. ','FE ','CHA ','CHB ','CHC ','CHD ','N A ','C1A ', + 'C2A ','C3A ','C4A ','CMA ','CAA ','CBA ','CGA ','O1A ', + 'O2A ','N B ','C1B ','C2B ','C3B ','C4B ','CMB ','CAB ', + 'CBB ','N C ','C1C ','C2C ','C3C ','C4C ','CMC ','CAC ', + 'CBC ','N D ','C1D ','C2D ','C3D ','C4D ','CMD ','CAD ', + 'CBD ','CGD ','O1D ','O2D ','1. ','CH3 ',15*'0 '/ C END C C C C ======================== SUBROUTINE ATMPOL(JVMAX) C ======================== C C---- establishes the vertices of the minimum required polyhedron C around the origin atom ao. C data taken from the ordered atom list in arrays ibox and box. C the results are loaded into the arrays ipolyh(i,j) i=1,3 (id C numbers of atoms in box contributing to the vertex j), and C polyhn(i,j) i=1,3 (xyz coordinates of vertex j). C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM,SAMSID,SAME INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) COMMON/FMRTST/NCALLS(25),JLEVEL,IMSG,IERROR(2,5),ICHK,JCHK DIMENSION P1(3),PLN(4) DIMENSION JELIM(20),LINES(2,60),ITEMP(4,100) EQUIVALENCE(TEMP(1,1),ITEMP(1,1)) save C C---- initialize error variables C JLEVEL=JLEVEL+1 C DO 90 J=JLEVEL,5 DO 90 I=1,2 90 IERROR(I,J)=0 C IERROR(1,JLEVEL)=17 C DO 91 I=1,100 91 VRTATM(I)=.FALSE. C C---- load ipolyh and polyhn with vertices from 4 hypothetical C atoms contained in box(i,j=97-100). C CALL VERTEX(97,98,99,1) CALL VERTEX(98,99,100,2) CALL VERTEX(99,100,97,3) CALL VERTEX(100,97,98,4) C C---- ipolyh(i,j) is loaded with the index numbers of atoms in box C whose planes intersect at vertex j. polyhn(i,j) is loaded C with the xyz coordinates of vertex j. C DO 100 I=97,100 100 VRTATM(I)=.TRUE. C C---- jvmax= total number of vertices included in polyhn C jatadd= number in box of atom whose plane is to produce C new vertices. C nelim= number of vertices in present list to be removed by C the plane of jatadd. C jelim(i)= index in ipolyh and polyhn of a vertex to be C eliminated from the list. i=1,nelim. C JVMAX=4 C C---- start checking vertices against rest of atoms in box. C 495 CONTINUE C C DO 500 I=1,20 500 JELIM(I)=0 C NELIM=0 JVERT=1 C C---- cycle thru all atoms in box or to next atom with plane that C eliminates 1 or more vertices. C DO 560 JJ=1,NASPHR C C---- skip atom if in vertex list or previously checked. C IF(VRTATM(JJ)) GO TO 560 C C---- load plane for this atom C DO 510 I=1,4 510 PLN(I)=BOX(I+1,JJ) C DD1 = 0.0 DD2 = 0.0 C C---- check each vertex against this plane C DO 550 J=1,JVMAX C C---- load vertex C DO 520 I=1,3 520 P1(I)=POLYHN(I,J) C IERROR(2,JLEVEL)=525 C C---- check vertex C 525 SAME=SAMSID(PLN,AO,P1,DD1,DD2) IF(SAME) GO TO 530 IF(JVERT.GT.20) GO TO 590 C C---- identify vertex to be eliminated C JELIM(JVERT)=J JVERT=JVERT+1 530 IF(J.NE.JVMAX) GO TO 550 IF(JELIM(1).EQ.0) GO TO 555 C C---- identify plane providing new vertex or vertices. C JATADD=JJ C C---- flag jj as a tested atom C VRTATM(JJ)=.TRUE. NELIM=JVERT-1 GO TO 600 550 CONTINUE C C---- flag jj as a tested atom. C 555 VRTATM(JJ)=.TRUE. 560 CONTINUE C C---- if nelim=0 there are no intersecting planes and no vertices C to be eliminated. C GO TO 800 590 IERROR(2,JLEVEL)=590 IMSG=8 CALL ERRORX GO TO 2000 C C---- remove and add vertices as indicated. C C---- load lines(2,30) with atom pairs derived from ipolyh(i,j) C at vertices j given by contents of jelim. C 600 DO 660 J=1,NELIM JE=JELIM(J) J3=3*J J2=J3-1 J1=J3-2 IF(IPOLYH(1,JE).GT.IPOLYH(2,JE)) GO TO 610 LINES(1,J1)=IPOLYH(1,JE) LINES(2,J1)=IPOLYH(2,JE) GO TO 620 610 LINES(1,J1)=IPOLYH(2,JE) LINES(2,J1)=IPOLYH(1,JE) 620 IF(IPOLYH(1,JE).GT.IPOLYH(3,JE)) GO TO 630 LINES(1,J2)=IPOLYH(1,JE) LINES(2,J2)=IPOLYH(3,JE) GO TO 640 630 LINES(1,J2)=IPOLYH(3,JE) LINES(2,J2)=IPOLYH(1,JE) 640 IF(IPOLYH(2,JE).GT.IPOLYH(3,JE)) GO TO 650 LINES(1,J3)=IPOLYH(2,JE) LINES(2,J3)=IPOLYH(3,JE) GO TO 660 650 LINES(1,J3)=IPOLYH(3,JE) LINES(2,J3)=IPOLYH(2,JE) 660 CONTINUE C C---- check each pair in lines(i=1-2,j) with rest of array ahead C of it. if any identity is found replace the pair jj with C zeros. otherwise compute new vertex. C C---- jvmax1 identifies the first empty slot in polyhn beyond block C of vertices now listed and before present elimination cycle. C JVMAX1=JVMAX+1 JE=1 NELIM3=3*NELIM C DO 710 JJ=1,NELIM3 C C---- has the pair jj been matched previously. C IF(LINES(1,JJ).EQ.0) GO TO 710 JJ1=JJ+1 IF(JJ1.GT.NELIM3) GO TO 680 DO 670 I=JJ1,NELIM3 IF(LINES(1,JJ).NE.LINES(1,I)) GO TO 670 IF(LINES(2,JJ).NE.LINES(2,I)) GO TO 670 LINES(1,I)=0 LINES(2,I)=0 GO TO 710 670 CONTINUE C C---- if the program gets thru the previous statement then pair jj C is unique. proceed to form new vertex with this line and C plane associated with jatadd. C 680 IF(JE.GT.NELIM) GO TO 690 JVERT=JELIM(JE) GO TO 700 690 IF(JVERT.GE.JVMAX1) GO TO 695 JVERT=JVMAX1 GO TO 700 695 JVERT=JVERT+1 700 JE=JE+1 JA=LINES(1,JJ) JB=LINES(2,JJ) IERROR(2,JLEVEL)=705 705 CALL VERTEX(JA,JB,JATADD,JVERT) IF(ICHK.EQ.-1) GO TO 2000 710 CONTINUE JVMAXC=JVMAX-NELIM+JE-1 C C---- all entries in ilines have now been checked and all vertices C removed or added as required. if number added is less than C number removed then some illegal entries in ipolyh and polyhn C need to be removed and the arrays contracted. C 720 IF(JE.GT.NELIM) GO TO 770 C C---- je is index in jelim(i) of next vertex to be eliminated. C JE1=JELIM(JE) JE=JE+1 IF(JE.GT.NELIM) GO TO 735 C DO 730 J=JE,NELIM JE2=JELIM(J) DO 730 I=1,3 IPOLYH(I,JE2)=0 730 POLYHN(I,JE2)=0.0 C C---- contract arrays ipolyh and polyhn. C 735 J1=JE1+1 DO 750 J=J1,JVMAX IF(IPOLYH(1,J).EQ.0) GO TO 750 DO 740 I=1,3 IPOLYH(I,JE1)=IPOLYH(I,J) 740 POLYHN(I,JE1)=POLYHN(I,J) JE1=JE1+1 750 CONTINUE C C---- clear high ends of contracted arrays. C DO 760 J=JE1,JVMAX DO 760 I=1,3 IPOLYH(I,J)=0 760 POLYHN(I,J)=0.0 C C---- reset jvmax to last entry in new array. C go back and check atom list in box C 770 JVMAX=JVMAXC GO TO 495 800 CONTINUE C C---- set vrtatm=.true. for all atoms which appear C in final ipolyh list. C DO 900 J=1,100 900 VRTATM(J)=.FALSE. DO 910 J=1,JVMAX DO 910 I=1,3 JJ=IPOLYH(I,J) 910 VRTATM(JJ)=.TRUE. C C---- check if any of vertices 97-100 are in final list . if yes C return error C DO 920 J=97,100 IF(.NOT.VRTATM(J)) GO TO 920 ICHK=-1 C C---- error write C WRITE(3,5) 5 FORMAT(' ONE OR MORE OF VERTICES 97-100 STILL IN POLYHN') C GO TO 2000 920 CONTINUE C 2000 JLEVEL=JLEVEL-1 C RETURN END C C C C ====================== SUBROUTINE BXLOAD(J,F) C ====================== C C---- loads constant factor, d,for equation for a plane C into box(5,j). a*x + b*y + c*z + d = 0. C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) C SAVE C SUMNP=0.0 C DO 100 K=1,3 K1=K+1 DF=AO(K) + F*BOX(K1,J) 100 SUMNP=SUMNP + BOX(K1,J)*DF C BOX(5,J)= -SUMNP RETURN END C C C C =============================== SUBROUTINE DET3C(A,B1,B2,B3,B4) C =============================== C C---- a is a 3x4 array. b1 is the value of the determinant C composed of columns 2,3,4; b2 of columns 1,3,4; b3 of columns C 1,2,4; b4 of columns 1,2,3. C DIMENSION A(3,4) C A12=A(2,1)*A(3,2)-A(2,2)*A(3,1) A13=A(2,1)*A(3,3)-A(2,3)*A(3,1) A14=A(2,1)*A(3,4)-A(2,4)*A(3,1) A23=A(2,2)*A(3,3)-A(2,3)*A(3,2) A24=A(2,2)*A(3,4)-A(2,4)*A(3,2) A34=A(2,3)*A(3,4)-A(2,4)*A(3,3) B1=A23*A(1,4)-A24*A(1,3)+A34*A(1,2) B2=A13*A(1,4)-A14*A(1,3)+A34*A(1,1) B3=A12*A(1,4)-A14*A(1,2)+A24*A(1,1) B4=A12*A(1,3)-A13*A(1,2)+A23*A(1,1) END C C C C =============================== SUBROUTINE DET3R(A,B1,B2,B3,B4) C =============================== C C---- a is a 4x3 array. b1 is the value of the determinant C composed of rows 2,3,4; b2 of rows 1,3,4; b3 of rows 1,2,4; C b4 of rows 1,2,3. C DIMENSION A(4,3) C A12=A(1,1)*A(2,2)-A(2,1)*A(1,2) A13=A(1,1)*A(3,2)-A(3,1)*A(1,2) A14=A(1,1)*A(4,2)-A(4,1)*A(1,2) A23=A(2,1)*A(3,2)-A(3,1)*A(2,2) A24=A(2,1)*A(4,2)-A(4,1)*A(2,2) A34=A(3,1)*A(4,2)-A(4,1)*A(3,2) B1=A23*A(4,3)-A24*A(3,3)+A34*A(2,3) B2=A13*A(4,3)-A14*A(3,3)+A34*A(1,3) B3=A12*A(4,3)-A14*A(2,3)+A24*A(1,3) B4=A12*A(3,3)-A13*A(2,3)+A23*A(1,3) END C C C ================== SUBROUTINE ERRORX C ================= C C---- this subroutine receives all error calls from all other C subroutines and the main program. the array ierror(i,j) C keeps track of the heirarchy of calling programs. j=1 C identifies the main program, j=2,--identify each subsequently C called subroutine in sequence with the last entry being the C routine calling error. i=1 contains a code number providing C the name of the subroutine required by j. i=2 contains the C statement number that produced the error in the calling pgm. C C imsg is a code number for the message to be printed. C jchk is disregarded at this time C where used ichk is set in program detecting and actual or C possible error and is tested at higher levels. C ichk=0= no error C ichk=1=error found but program to continue with or without C message C ichk=-1= error found index used to abort program if necessary C after return from error. C C the only subroutines with /fmrtst/ at this time are: C samsid,inside,point,plane,vertex,main,error,atvol C only point,plane,atvol call error directly C COMMON/FMRTST/NCALLS(25),JLEVEL,IMSG,IERROR(2,5),ICHK,JCHK CHARACTER *4 NAME(50) LOGICAL VERBOSE COMMON /LOGCOM/ VERBOSE SAVE C DATA NAME/' ',' ','DET3','C ','DET3','R ', 1 'RELR','AD ','SAMS','ID ','CENT','RD ', 2 'DPTP','LN ','INSI','DE ','POIN','T ', 3 'PLAN','E ','VOLT','ET ','SELE','CT ', 4 'REST','AK ','ORDE','R ','VERT','EX ', 5 'NCAL','LS ','ERRO','R ','ATVO','L ', 6 'PRIN','T1 ','MAIN',' ',10*' '/ C WRITE(3,1) 1 FORMAT(/,' ERRORX ROUTINE CALLED',/) C IF (VERBOSE) THEN DO 100 J=1,5 INAME1=(IERROR(1,J)*2)+1 INAME2=INAME1+1 WRITE(3,2) NAME(INAME1),NAME(INAME2),IERROR(2,J) 2 FORMAT(1X,5(A4,A4,2X,I5),/) 100 CONTINUE END IF 200 CONTINUE C GO TO (1000,1000,300,400,500,600,700,800),IMSG C 300 WRITE(3,3) 3 FORMAT(' THE PLANE IS INDETERMINATE. THE THREE DEFINING POINTS ARE 1 COLINEAR.') GO TO 1000 C 400 WRITE(3,4) 4 FORMAT(' NO POINT OF INTERSECTION. PLANES ARE PARALLEL, INTERSECT' 1,/,' IN A LINE OR INTERSECT IN PAIRS TO GIVE 3 LINES.') GO TO 1000 C 500 GO TO 1000 C 600 WRITE(3,6) 6 FORMAT(' THE PROGRAM CANNOT FIND 3 ATOMS IN BOX TO DEFINE A ',/, 1' PLANE MORE THAN 0.5 A FROM AO, AND ON OPPOSITE SIDE OF AO',/, 2' FROM ATOM 1 IN BOX. PICK NEXT ATOM IN INPUT LIST.') GO TO 1000 C 700 WRITE(3,7) 7 FORMAT(' THE PROGRAM CANNOT FIND 3 ATOMS IN BOX TO COMBINE',/, 1' WITH ATOM NUMBER 1 TO DEFINE A TETRAHEDRON WHICH SURROUNDS',/, 2' THE ORIGIN ATOM AO. PICK NEXT ATOM IN INPUT LIST.') GO TO 1000 C 800 WRITE(3,8) 8 FORMAT(' THE NUMBER OF VERTICES IN THE PRESENT POLYHEDRON TO ', 1/,' BE ELIMINATED BY THE PLANE ASSOCIATED WITH THE NEXT ATOM ', 2/,' IN BOX IS GREATER THAN 10, THE PRESENT MAXIMUM OF NELIM.',/, 3' PICK NEXT ATOM IN INPUT LIST.') C 1000 RETURN END C C C C ================== SUBROUTINE FILBOX C ================== C C---- loads arrays box and ibox with coordinates and ids of atoms C within the sphere of nominal radius arange C centered on the position of the origin atom ao. C C temp(1,kk)= box(1,j)= distance of atom (kk or j) from origin atom ao. C temp(2,kk)= box(2,j)= delta x ( kk or j from ao ) C temp(3,kk)= box(3,j)= delta y ( kk or j from ao ) C temp(4,kk)= box(4,j)= delta z ( kk or j from ao ) C box(5,j)= constant factor in equation of plane (d). C C itemp(1,kk)= ibox(1,j)= serial number of atom or 'solvent' in coord lst C itemp(2,kk)= ibox(2,j)= 0 for protein atom, 1 for solvent 'atom' C C equation for plane: a*x + b*y + c*z + d = 0 C C a,b,c are the direction components of the plane normal C x,y,z are the coordinates of any point in the plane C but in this subroutine it is the point of intersection of the C plane and the interatomic vector to which it is normal. C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2, + AOSEQ,AONUM,BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100), + VRTATM(100),VERVOL(100), 1 BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM DIMENSION AL(3),ITEMP(2,100) DIMENSION NX(3),XC(3),XCL(3),XCH(3) C SAVE C C---- default van der waals radius C VDWDFT=1.70 C NX(1)=NINT ((AO(1)-XMIN)/CE ) NX(2)=NINT ((AO(2)-YMIN)/CE ) NX(3)=NINT ((AO(3)-ZMIN)/CE ) C XC(1)=(FLOAT(NX(1))+0.5)*CE+XMIN XC(2)=(FLOAT(NX(2))+0.5)*CE+YMIN XC(3)=(FLOAT(NX(3))+0.5)*CE+ZMIN C DO 30 J=1,3 XCL(J)=XC(J)-CE-0.1 30 XCH(J)=XC(J)+CE+0.1 C C nx(i) = 3 less than actual indices of icl cube containing ao C in shell. C xc(i) = floating point coordinates of this cube. C xcl(i) and xch(i) define the limits of the 26 surrounding cubes. C the cube edge is ce, the same number used in defining the C original solvent lattice. C DO 888 JDO=1,100 ITEMP(1,JDO) = 0 ITEMP(2,JDO) = 0 888 CONTINUE C RNGSQ1=ARANGE*ARANGE NTEST=NPAT+NSOL KK=1 I=0 50 I=I+1 IF(I.GT.NTEST) GO TO 305 C C---- delete all protein atoms with key=-1 C IF(I.LE.NPAT.AND.KEY(I).LT.0) GO TO 300 AL(1)=X(I)-AO(1) AL(2)=Y(I)-AO(2) AL(3)=Z(I)-AO(3) SUMSQ=0.0 DO 100 J=1,3 100 SUMSQ=SUMSQ+AL(J)*AL(J) JSEQN=SEQNUM(I) C C---- lstres is used as the maximum sequence number of protein C atoms in coordinate list. higher numbers refer to C hypothetical solvent molecules C IF(JSEQN.GT.LSTRES) GO TO 150 IF(SUMSQ.GT.RNGSQ1) GO TO 300 C C---- eliminate ao itself from list C IF(SUMSQ.LT.0.0001) GOTO 300 GO TO 250 150 CONTINUE C C---- load solvent atom in temp if it is any of 26 nearest C neighbors of spherical ao. C IF(X(I).LT.XCL(1).OR.X(I).GT.XCH(1)) GO TO 300 IF(Y(I).LT.XCL(2).OR.Y(I).GT.XCH(2)) GO TO 300 IF(Z(I).LT.XCL(3).OR.Z(I).GT.XCH(3)) GO TO 300 C 250 CONTINUE C C---- see notes in header to subroutine C TEMP(1,KK)=SQRT(SUMSQ) TEMP(2,KK)=AL(1) TEMP(3,KK)=AL(2) TEMP(4,KK)=AL(3) ITEMP(1,KK)=I IF(JSEQN.GT.LSTRES) ITEMP(2,KK)=1 KK=KK+1 300 CONTINUE GO TO 50 305 CONTINUE KK=KK-1 NASPHR=KK C C---- total number of atoms in sphere around ao = nasphr.(symbol C kk is used internally in this subroutine.) C VDIAM = (AOVDW + RDELTA) * 2.0 AOVDW2=AOVDW**2 C C---- if second atom is solvent molecule, arbitrarily set C dividing plane at the van der waals surface of ao. C or at a surface rdelta outside of the actual vander C waals surface. C DO 309 I = 1,KK IF (ITEMP (2,I) .EQ. 0) GO TO 309 RATIO = VDIAM/TEMP (1,I) TEMP (1,I) = VDIAM DO 308 K = 2,4 308 TEMP (K,I) = TEMP (K,I) * RATIO 309 CONTINUE C C---- put atoms in box in order of increasing distance from ao. C C ll = index for array box C l = index for array temp C LL=1 L=KK 310 CONTINUE DNEXT=1000.0 DO 320 I=1,L IF(DNEXT.LE.TEMP(1,I)) GO TO 320 DNEXT=TEMP(1,I) INEXT=I 320 CONTINUE DO 340 J=1,4 340 BOX(J,LL)=TEMP(J,INEXT) IBOX(1,LL)=ITEMP(1,INEXT) IBOX(2,LL)=ITEMP(2,INEXT) 350 LL=LL+1 C C---- remove selected atom from temp. C IF(INEXT.EQ.L) GO TO 380 DO 360 J=1,4 360 TEMP(J,INEXT)=TEMP(J,L) ITEMP(1,INEXT)=ITEMP(1,L) ITEMP(2,INEXT)=ITEMP(2,L) 380 IF(L.EQ.1) GO TO 390 L=L-1 GO TO 310 390 CONTINUE C GO TO (400,500,600), METHOD C C---- load box(5,j) by normal voronoi procedure C ( richards method a ) (i.e. dividing plane bisects C interatomic vectors) C 400 F=0.5 DO 440 J=1,KK 440 CALL BXLOAD(J,F) GO TO 700 C C---- load box(5,j) by richards method b (i.e. dividing C plane positioned by radii ratios) C 500 DO 540 J=1,KK JID=IBOX(1,J) C C---- if j is a 'solvent' atom, bypass covalent check C IF(SEQNUM(JID).GT.LSTRES) THEN F=0.5 GO TO 530 ENDIF C C---- determine if j is covalently bonded to ao or not. C SEQDIF=SEQNUM(JID) - AOSEQ C C---- get residue type number and atom number C CALL GETNUMX (JID,NRES,JNUM) C C---- error write C IF(NRES.EQ.0.OR.JNUM.EQ.0) THEN WRITE(3,5000) NRES,JNUM,J 5000 FORMAT(' FILBOX: NRES =',I5,' NATM =',I5,' BOX ATOM =',I5, + ' DEFAULT VDW RADIUS USED') GO TO 515 END IF C IF(AONUM.EQ.2.AND.SEQDIF.EQ.-1.AND.JNUM.EQ.3) GO TO 520 IF(AONUM.EQ.3.AND.SEQDIF.EQ. 1.AND.JNUM.EQ.2) GO TO 520 C C---- set up bond code C BNDCOD=AONUM*100 + JNUM IF(AONUM.GT.JNUM) BNDCOD=JNUM*100 + AONUM C C---- check bond code list C IN=INDX2(NRES) INCR=IBND1(IN) IB=IN+1 IE=IN+INCR DO 510 I=IB,IE IF(BNDCOD.NE.IBND1(I)) GO TO 510 GO TO 520 510 CONTINUE C C---- set up plane based on division by vanderwaals radii. C F=(BOX(1,J) + AOVDW - RVDW(JID))/(2.*BOX(1,J)) CALL BXLOAD(J,F) GO TO 540 515 F=(BOX(1,J) + AOVDW - VDWDFT)/(2.*BOX(1,J)) CALL BXLOAD(J,F) GO TO 540 C C---- set up plane based on division by covalent radii. C 520 F= AOCOV/(AOCOV+RCOV(JID)) 530 CALL BXLOAD(J,F) 540 CONTINUE GO TO 700 C C---- load box(5,j) by the radical plane procedure C (see gellatly and finney j.mol.biol.161,305-322 (1982)) C 600 DO 640 J=1,KK IF(IBOX(2,J).EQ.0) GO TO 620 C C---- when atom j is a 'solvent'atom, the plane is drawn C tangent to the vdw sphere of ao ( see earlier statements C in this subroutine which reset temp((1-4,j) for this C eventuality.) C F=0.5 CALL BXLOAD(J,F) GO TO 640 C C---- when atom j is a protein atom, the radical plane procedure C provides the locus of points equidistant from the tangents to C the vanderwaals spheres. C 620 R2SQ=RVDW(IBOX(1,J))**2 DSQ=BOX(1,J)**2 F=(DSQ + AOVDW2 - R2SQ)/(2*DSQ) CALL BXLOAD(J,F) 640 CONTINUE C C---- load box positions 97-100 with hypothetical tetrahedron C surrounding ao .to be used as starting polyhn in atpol. C 700 CONTINUE BOX(2,97) =+10.0 BOX(2,98) =+10.0 BOX(2,99) =-10.0 BOX(2,100)=-10.0 BOX(3,97) =+10.0 BOX(3,98) =-10.0 BOX(3,99) =+10.0 BOX(3,100)=-10.0 BOX(4,97) =+10.0 BOX(4,98) =-10.0 BOX(4,99) =-10.0 BOX(4,100)=+10.0 C DO 800 J=97,100 800 BOX(1,J)=17.3205 C F=0.5 DO 820 J=97,100 820 CALL BXLOAD(J,F) RETURN END C C C C =============================== SUBROUTINE GETNUMX(J,NRES,NATM) C =============================== C C---- returns residue type number,nres, and atom type number C within that residue, natm, for an atom identified by C the serial number j. C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) C CHARACTER ATMTMP*4,RESTMP*4 C SAVE NRES=0 NATM=0 RESTMP=RES3(J) DO 100 I=1,30 IF(RESTMP.NE.RESTYP(I)) GO TO 100 NRES=I GO TO 150 100 CONTINUE C WRITE(3,1000) RESTMP 1000 FORMAT(/,' ERROR: RESIDUE NOT FOUND IN GETNUMX:',2X,A4) C 150 ATMTMP=ATM(J) IF(ATMTMP.EQ.'CA ') THEN NATM=1 GO TO 500 ELSE IF(ATMTMP.EQ.'N ') THEN NATM=2 GO TO 500 ELSE IF(ATMTMP.EQ.'C ') THEN NATM=3 GO TO 500 ELSE IF(ATMTMP.EQ.'O ') THEN NATM=4 GO TO 500 ENDIF C IN=INDX1(NRES) INCR=KATNAMA(IN) IB=IN+1 IE=IN+INCR DO 300 I=IB,IE IF(ATMTMP.NE.ATNAMA(I)) GO TO 300 NATM=I-IB+1 IF(NRES.LT.22) NATM=NATM+4 GO TO 500 300 CONTINUE C WRITE(3,1100) J,ATMTMP 1100 FORMAT(/,' ERROR: ATM NOT FOUND IN GETNUMX. SERIAL # =',I5, 1 ' NAME =',2X,A4) C 500 RETURN C END C C C C =========================== SUBROUTINE POINT (PLANES,P) C =========================== C C---- planes is a 4x3 array representing equations of 3 planes. C the subroutine calculates their point of intersection p. C planes(1-3,j)= direction components of plane j, C plane(4,j)= constant term of plane j. C the value of ichk = 0 if the planes have a point of C intersection, ichk=-1 if no point of intersection (i.e. C planes are parallel, intersect in a line, or intersect in C 3 lines in pairs) C COMMON/FMRTST/NCALLS(25),JLEVEL,IMSG,IERROR(2,5),ICHK,JCHK DIMENSION PLANES(4,3), P(3) LOGICAL VERBOSE COMMON /LOGCOM/ VERBOSE SAVE C C---- initialize error variables C JLEVEL=JLEVEL+1 DO 90 J=JLEVEL,5 DO 90 I=1,2 90 IERROR(I,J)=0 IERROR(1,JLEVEL)=8 ICHK=0 C CALL DET3R(PLANES,B1,B2,B3,B4) IF(B4.NE.0) GO TO 100 C IMSG=4 IERROR(2,JLEVEL)=50 50 CALL ERRORX ICHK=-1 IF (VERBOSE) THEN WRITE(3,2) 2 FORMAT(6X,'NX',12X,'NY',12X,'NZ',12X,'D',/) DO 60 J=1,3 WRITE(3,3)(PLANES(I,J),I=1,4) 3 FORMAT(4F12.6) 60 CONTINUE END IF GO TO 1000 100 P(1)=-B1/B4 P(2)=B2/B4 P(3)=-B3/B4 1000 JLEVEL=JLEVEL-1 RETURN END C C C C ======================== SUBROUTINE POLVOL(JVMAX) C ======================== C C---- calculates the volume of the irregular polyhedron whose C vertices are listed in the arrays ipolyh and polyhn. C C---- the procedure divides each face of the polyhedron into C triangles. computes the volume of the little tetrahedron C formed by each triangle and the central atom ao. then sums C all the tetrahedral volumes to get the volume of the C polyhedron. C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) C DIMENSION VERT(3,4),ITEMP(4,100) C EQUIVALENCE (TEMP(1,1),ITEMP(1,1)) C SAVE C VOLUME=0.0 C C----- load the common vertex ao into array vert(i=1-3,1) C DO 100 I=1,3 100 VERT(I,1)=AO(I) C DO 1000 NXTAT=1,100 C C---- check that the atom to be used is in the vertex list. C IF(.NOT.VRTATM(NXTAT)) GO TO 1000 C C---- load itemp with all vertices from ipolyh containing nxtat. C JJ=1 DO 200 J=1,JVMAX DO 180 I=1,3 IF(IPOLYH(I,J).NE.NXTAT) GO TO 180 ITEMP(1,JJ)=IPOLYH(1,J) ITEMP(2,JJ)=IPOLYH(2,J) ITEMP(3,JJ)=IPOLYH(3,J) ITEMP(4,JJ)=J JJ=JJ+1 GO TO 200 180 CONTINUE 200 CONTINUE JN=JJ-1 IF(JN.EQ.0) GO TO 1500 C C---- all entries in itemp(i,jj) contain nxtat in locations 1,2, or C 3 for each value of jj. C jn=total number of vertices loaded into itemp. C C---- first entry in itemp is the second vertex for this round C of calculation C ID2=ITEMP(4,1) DO 210 I=1,3 210 VERT(I,2)=POLYHN(I,ID2) C C----- determine either of the other two atoms,=n2, in first entry C that is not nxtat. C DO 220 I=1,3 N2=ITEMP(I,1) IF(N2.NE.NXTAT) GO TO 240 220 CONTINUE C C---- identify third vertex using n2 C 240 DO 250 J=2,JN J3= J DO 250 I=1,3 IF(ITEMP(I,J).EQ.N2) GO TO 260 250 CONTINUE C C---- load third vertex into vert. C 260 ID3=ITEMP(4,J3) DO 270 I=1,3 270 VERT(I,3)=POLYHN(I,ID3) NLAST=N2 JNOW=J3 IDNOW=ID3 IX=1 C C---- identify atom which is not nxtat or nlast in vertex jnow. C 300 DO 380 I=1,3 NNOW=ITEMP(I,JNOW) IF(NNOW.EQ.NXTAT) GO TO 380 IF(NNOW.EQ.NLAST) GO TO 380 GO TO 400 380 CONTINUE C C---- Identify vertex jnext using nnow. C 400 DO 450 J=2,JN JNEXT=J DO 450 I=1,3 IF(ITEMP(I,J).NE.NNOW) GO TO 450 IF(ITEMP(4,J).EQ.IDNOW) GO TO 450 GO TO 460 450 CONTINUE C C---- load vertex jnext into array vert. C 460 IDNEXT=ITEMP(4,JNEXT) JALT=3+MOD(IX,2) DO 470 I=1,3 470 VERT(I,JALT)=POLYHN(I,IDNEXT) C C---- calculate volume of this tetrahedral segment. C 660 VOLI=VOLTET(VERT) VOLUME=VOLUME+VOLI IF((IX+2).EQ.JN) GO TO 1000 NLAST=NNOW JNOW=JNEXT IDNOW=IDNEXT IX=IX+1 GO TO 300 1000 CONTINUE GO TO 2000 C C---- error write C 1500 WRITE(3,6) 6 FORMAT(' ERROR IN VRTATM OR IPOLYH SUPPLIED TO POLVOL') 2000 RETURN END C C C C ===================== SUBROUTINE PRNTBX(ID) C ===================== C C---- prints out current contents of arange,id,ao(3), C ibox(jj),box(i=1-11,jj) where jj=1-nacube or 1-nasphr. C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) LOGICAL VERBOSE COMMON /LOGCOM/ VERBOSE C SAVE C WRITE(3,1) ARANGE,ID,AO(1),AO(2),AO(3) 1 FORMAT( //,' MAXIMUM RADIUS=',F8.3,', ID=', I5,' AO: X=',F8.3, +', Y=',F8.3,', Z=',F8.3) WRITE(3,2) 2 FORMAT(/,' IBOX(1-2,J),BOX(I,J),I=1,5,J=1-NASPHR') C IF (VERBOSE) THEN WRITE(3,3) 3 FORMAT(/,' NO.',5X,'ID',3X,' SOLV',5X, + 'DIST',7X,'NX', + 8X,'NY',8X,'NZ',8X,'D',/) DO 100 J=1,NASPHR 100 WRITE(3,4)J,IBOX(1,J),IBOX(2,J),(BOX(I,J),I=1,5) 4 FORMAT(I4,I10,I5,5F10.5) WRITE(3,5) 5 FORMAT (//) DO 200 J = 97,100 200 WRITE(3,4)J,IBOX(1,J),IBOX(2,J),(BOX(I,J),I=1,5) WRITE (3,5) END IF C END C C C C ======================== SUBROUTINE PRNTPN(JVMAX) C ======================== C C---- prints out current contents of arrays ipolyh and polyhn. C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) LOGICAL VERBOSE COMMON /LOGCOM/ VERBOSE C SAVE C WRITE(3,5) 5 FORMAT(//,' IPOLYH(I,J),POLYHN(I,J),I=1-3,J=1-JVMAX') WRITE(3,6) 6 FORMAT(/,' NO.',4X,'JA',8X,'JB',8X,'JC',8X,'X',9X,'Y',9X,'Z',/) C IF (VERBOSE) THEN DO 200 J=1,JVMAX 200 WRITE(3,7)J,(IPOLYH(I,J),I=1,3),(POLYHN(I,J),I=1,3) 7 FORMAT(I4,3(4X,I2,4X),3F10.5) END IF C END C C C C ============================================= LOGICAL FUNCTION SAMSID(PLANE1,P1,P2,DD1,DD2) C ============================================= C C---- determines whether points p1 and p2 are on same side of a C given plane,'plane1'. yes=.true. , no=.false. C plane1(1-3)= direction components of plane normal C plane1(4)= constant term C p1(1-3),p2(1-3) are the xyz coordinates of the two points C C---- note that if one of the points tested by samsid is in the C plane, the result is indeterminate.the function is set C equal to .true. and ichk is set equal to 1. normally ichk=0 C COMMON/FMRTST/NCALLS(25),JLEVEL,IMSG,IERROR(2,5),ICHK,JCHK DIMENSION PLANE1(4),P1(3),P2(3) SAVE C ICHK=0 D1=0.0 D2=0.0 C C---- calculate constant term for planes parallel to plane but C passing through p1 and p2 C DO 100 I = 1,3 100 D1 = D1-PLANE1 (I) * P1(I) DD1 = PLANE1(4) -D1 140 DO 150 I = 1,3 150 D2 = D2-PLANE1(I)*P2(I) DD2 = PLANE1(4)-D2 C C---- test constants against that of plane C DD1SQ = DD1 * DD1 DD2SQ = DD2 * DD2 IF (DD1SQ .LT. 0.0001) GO TO 300 IF (DD2SQ .LT. 0.0001) GO TO 300 IF(DD1*DD2) 200,300,400 200 SAMSID=.FALSE. RETURN C 300 ICHK=1 400 SAMSID=.TRUE. RETURN END C C C C ================ SUBROUTINE SHELL C ================ C C---- program to put a rough solvent shell around a protein molecule. C the shell is based on a cubic lattice of edge dimension ce. C the shell, in the format of the protein coordinate files, C is output as a separate file on lun 4. this file, either old C or new, is read in by volume for use in the subsequent C calculations. C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) C DIMENSION NCOUNT(10),NC(30) CHARACTER*2 NC CHARACTER*4 XRES3,XATM C SAVE C DO 50 I=1,MAXSIZ 50 ICL(I)= 0 C C---- the array icl is set up as a one dimensional array but C is thought of as a three dimensional array in describing C the operation of the program. C set up the maximum dimensions of the required array in C the three axial directions. C IMAX=NINT((XMAX-XMIN)/CE+5.) JMAX=NINT((YMAX-YMIN)/CE+5.) KMAX=NINT((ZMAX-ZMIN)/CE+5.) IM=IMAX JIM=JMAX*IMAX C IF(IMAX*JMAX*KMAX.GT.MAXSIZ) GO TO 5000 C C---- xm2,ym2,zm2 are the center of the first cube in the icl array C XM2=XMIN+CE2 YM2=YMIN+CE2 ZM2=ZMIN+CE2 C C---- load atoms into cubes in array icl(l) .smallest starting C index for each axis =3. C DO 150 I = 1,NPAT IF(KEY(I).LT.0) GO TO 150 IN=NINT((X(I)-XMIN)/CE+3.) JN=NINT((Y(I)-YMIN)/CE+3.) KN=NINT((Z(I)-ZMIN)/CE+3.) L= (KN-1)*JIM + (JN-1)*IM + IN ICL(L)=ICL(L)+1 150 CONTINUE C C---- icl entries now equal zero or the number of protein atoms C in that cube. C LCNT = 0 DO 200 K = 1,KMAX KM = (K-1)*JIM DO 200 J = 1,JMAX JM = (J-1)*IM DO 200 I = 1,IMAX L = KM + JM + I IF (ICL(L).LE.0) GO TO 200 LDEL=ICL(L) L1 = LCNT + 1 ICL(L) = L1 JSERN(L1) = LDEL + 10000*LDEL LCNT = L1 + LDEL 200 CONTINUE C C---- icl now contains pointer to the jsern location for cube(i,j,k) which C gives the total number of entries for that cube. atom id C numbers for that cube are to follow as next entries in jsern. C DO 300 I = 1,NPAT IF(KEY(I).LT.0) GO TO 300 IN=NINT((X(I)-XMIN)/CE+3.) JN=NINT((Y(I)-YMIN)/CE+3.) KN=NINT((Z(I)-ZMIN)/CE+3.) L = (KN-1)*JIM + (JN-1)*IM + IN INDX = ICL (L) J = INDX + MOD(JSERN(INDX),10000) JSERN (J) = I JSERN(INDX) = JSERN(INDX) -1 IF(MOD(JSERN(INDX),10000).GT.1000) GO TO 320 300 CONTINUE GO TO 350 C C---- error write C 320 CALL CCPERR(1, ' stop shell error at statement 320') C C---- reset jsern indx points C 350 CONTINUE DO 400 K = 1,KMAX KM = (K-1)*JIM DO 400 J = 1,JMAX JM = (J-1)*IM DO 400 I = 1,IMAX L = KM + JM + I IF (ICL(L).LE.0) GO TO 400 INDX = ICL(L) JSERN(INDX)=JSERN(INDX)/10000 400 CONTINUE C C----- load -1 into icl(l) if all neighbors are 0 or -1 C DO 700 K=1,KMAX KB = K-1 IF(KB.LT.1) KB=1 KE = K+1 IF(KE.GT.KMAX) KE=KMAX DO 700 J=1,JMAX JB = J-1 IF(JB.LT.1) JB=1 JE = J+1 IF(JE.GT.JMAX) JE=JMAX DO 700 I=1,IMAX IB = I-1 IF(IB.LT.1) IB=1 IE = I+1 IF(IE.GT.IMAX) IE=IMAX M=0 DO 600 KK=KB,KE KN=(KK-1)*JIM DO 600 JJ=JB,JE JN=(JJ-1)*IM DO 600 II=IB,IE LL = KN + JN + II IF(ICL(LL).LE.0) GO TO 600 M=M+1 600 CONTINUE IF(M.EQ.0) GO TO 650 GO TO 700 650 L= (K-1)*JIM + (J-1)*IM + I ICL(L)=-1 700 CONTINUE C C---- potential solvent molecules now restricted to positions C labelled 0. search for channel to surface in each case. C 841 IMAX1=IMAX-1 JMAX1=JMAX-1 KMAX1=KMAX-1 DO 1000 K=2,KMAX1 KM=(K-1)*JIM DO 1000 J=2,JMAX1 JM=(J-1)*IM DO 1000 I=2,IMAX1 L = KM + JM + I IF(ICL(L).NE.0) GO TO 1000 C C---- examine 26 sites around position i, j, k. C DO 950 KK=1,3 KN=K+KK-2 DO 950 JJ=1,3 JN=J+JJ-2 DO 950 II=1,3 IN=I+II-2 IF(II.EQ.2.AND.JJ.EQ.2.AND.KK.EQ.2) GO TO 950 LL=(KN-1)*JIM+(JN-1)*IM+IN IF(ICL(LL).LE.0) GO TO 950 C C---- test if any opposing pairs of lattice sites around target C atom are both occupied. if yes, label cube as a potential C cavity (-2). C INOP=IN+2*(2-II) JNOP=JN+2*(2-JJ) KNOP=KN+2*(2-KK) LL= (KNOP-1)*JIM + (JNOP-1)*IM + INOP IF(ICL(LL).LE.0) GO TO 950 ICL(L)=-2 GO TO 1000 950 CONTINUE 1000 CONTINUE C C---- check if any solvent or -2 positions are inside protein C vdw envelopes.if yes,set equal to -3 or -4. C DO 2200 I=1,NPAT IF(KEY(I).LT.0) GO TO 2200 R=RVDW(I) RSQ=R**2 IN=NINT((X(I)-XMIN)/CE+3.) JN=NINT((Y(I)-YMIN)/CE+3.) KN=NINT((Z(I)-ZMIN)/CE+3.) IS=IN-1 JS=JN-1 KS=KN-1 IE=IS+2 JE=JS+2 KE=KS+2 DO 2120 KK=KS,KE KM = (KK-1)*JIM DO 2120 JJ=JS,JE JM = (JJ-1)*IM DO 2120 II=IS,IE AX=XM2+CE*FLOAT(II-3) AY=YM2+CE*FLOAT(JJ-3) AZ=ZM2+CE*FLOAT(KK-3) DIFX=ABS(AX-X(I)) DIFY=ABS(AY-Y(I)) DIFZ=ABS(AZ-Z(I)) IF(DIFX.GT.R.OR.DIFY.GT.R.OR.DIFZ.GT.R) GO TO 2120 CHKSQ=DIFX**2+DIFY**2+DIFZ**2 IF(CHKSQ.GT.RSQ) GO TO 2120 LL = KM + JM + II IF(ICL(LL).EQ.0) ICL(LL)=-3 IF(ICL(LL).EQ.-2) ICL(LL)=-4 2120 CONTINUE 2200 CONTINUE C C---- set -2 positions that are in contact with solvent to -5. C residual -2 positions are thus truly interior cavities. C DO 2500 K=2,KMAX1 KM=(K-1)*JIM DO 2500 J=2,JMAX1 JM=(J-1)*IM DO 2500 I=2,IMAX1 L = KM + JM + I IF(ICL(L).NE.-2) GO TO 2500 C C---- examine 26 sites around position i, j, k. C DO 2400 KK=1,3 KN=K+KK-2 DO 2400 JJ=1,3 JN=J+JJ-2 DO 2400 II=1,3 IN=I+II-2 IF(II.EQ.2.AND.JJ.EQ.2.AND.KK.EQ.2) GO TO 2400 LL=(KN-1)*JIM+(JN-1)*IM+IN IF(ICL(LL).EQ.0) THEN ICL(L)=-5 GOTO 2500 ENDIF 2400 CONTINUE 2500 CONTINUE 2000 NZERO= 0 MIN2 = 0 MIN3 = 0 MIN4 = 0 MIN5 = 0 NPROT= 0 DO 1090 N=1,10 1090 NCOUNT(N)=0 DO 1100 K=1,KMAX KM = (K-1)*JIM DO 1100 J=1,JMAX JM = (J-1)*IM DO 1100 I=1,IMAX L = KM + JM + I JCL = ICL(L) IF(JCL.EQ.0) NZERO = NZERO+1 IF(JCL.EQ.-2) MIN2 = MIN2 + 1 IF(JCL.EQ.-3) MIN3 = MIN3 + 1 IF(JCL.EQ.-4) MIN4 = MIN4+1 IF(JCL.EQ.-5) MIN5 = MIN5 + 1 IF(JCL.GT.0) NPROT = NPROT + 1 IF(JCL.LE.0) GO TO 1100 N=JSERN(JCL) NCOUNT(N)=NCOUNT(N) + 1 1100 CONTINUE C C---- write to data output tape C WRITE(2,8) NZERO,MIN3,MIN2,MIN4,MIN5,NPROT,(NCOUNT(N),N=1,10) 8 FORMAT(/,' NUMBER OF SOLVENT MOLECULES REMAINING ( 0)=',I10,/, 1' NUMBER OF SOLVENT MOLECULES REMOVED (-3)=',I10/, 2' NUMBER OF INTERIOR CAVITIES REMAINING (-2)=',I10/, 3' NUMBER OF INTERIOR CAVITIES REMOVED (-4)=',I10/, 4' NUMBER OF SURFACE CAVITIES REMAINING (-5)=',I10/, 5' NUMBER OF CELLS WITH PROTEIN ATOMS (**)=',I10,/, 6' NUMBER OF CELLS WITH 1,2,3,ETC PROTEIN ATOMS = ',10I5,/) C C---- load solvent 'atoms' into output file on lun 4 C MAXSOL = LSTRES + 1 + (NZERO/100) NTOT = NPAT + NZERO + MIN2 + MIN5 WRITE(4,5100) MAXSOL,NZERO,MIN2,MIN5,NTOT 5100 FORMAT(4I5) MSOL = 0 IKEY = 1 XRES3 = 'HOH ' XATM = 'HOH ' DO 3000 K = 1, KMAX KM = (K-1)*JIM DO 3000 J = 1, JMAX JM = (J-1)*IM DO 3000 I = 1, IMAX L = KM + JM + I IF (ICL(L) .NE. 0) GO TO 3000 IRSDEL = MSOL/100 IRES = LSTRES + 1 + IRSDEL MSOL = MSOL + 1 LTST = NPAT + MSOL IF (LTST .GT. MAXATM) GO TO 4500 XX=XM2+CE*FLOAT(I-3) XY=YM2+CE*FLOAT(J-3) XZ=ZM2+CE*FLOAT(K-3) C C---- nchn,nfil,xvdw,xcov,xarea,yarea,zarea, are all dummy C variables set = 0 to make file record format C equivalent to input file format. C NCHN = 0 NFIL = 0 XVDW = 0 XCOV = 0 XAREA = 0 YAREA = 0 ZAREA = 0 WRITE(4,5200) IKEY,NCHN,NFIL,XATM,XRES3,IRES,XX,XY,XZ,XVDW, 1 XCOV,XAREA,YAREA,ZAREA 5200 FORMAT(1X,I2,7X,I2,1X,I2,1X,2A4,1X,I3,3F8.3, 1 2F5.2,2F6.1,F5.2) 3000 CONTINUE C C---- load cubes labeled -2 or -5 at end of output file C NMIN=0 IKEY = 1 XRES3 = 'DUM ' XATM = 'DUM ' DO 3100 K = 1, KMAX KM = (K-1)*JIM DO 3100 J = 1, JMAX JM = (J-1)*IM DO 3100 I = 1, IMAX L = KM + JM + I IF (ICL(L) .NE. -2.AND.ICL(L).NE.-5) GO TO 3100 IRSDEL = NMIN/100 IRES = MAXSOL + 1 + IRSDEL NMIN = NMIN + 1 LTST = NPAT + MSOL + NMIN C IF (LTST .GT. MAXATM) GO TO 4500 C XX=XM2+CE*FLOAT(I-3) XY=YM2+CE*FLOAT(J-3) XZ=ZM2+CE*FLOAT(K-3) WRITE(4,5200) IKEY,NCHN,NFIL,XATM,XRES3,IRES,XX,XY,XZ,XVDW, 1 XCOV,XAREA,YAREA,ZAREA 3100 CONTINUE GO TO 7000 C C---- error write C 4500 WRITE(3,4510) LTST 4510 FORMAT ('NUMBER OF PROTEIN + SOLVENT ATOMS GREATER THAN',I10) GO TO 6000 C C--- error write C 5000 CONTINUE ICLX = 0 ICLY = 0 ICLZ = 0 WRITE(6,5010) IMAX,JMAX,KMAX,ICLX,ICLY,ICLZ 5010 FORMAT(' ICL ARRAY TOO SMALL',/,' IMAX,JMAX,KMAX =',8I5,/, 1' ICLX,ICLY,ICLZ = ',3I5) 6000 STOP 7000 CONTINUE C C---- to get an image of icl compile 'shell' with d_lines. C the image is added on at the end of the shell output C file on lun 4. note that -3 and -4 positions are shown C as part of the protein (**). C DO 7020 K=1,KMAX WRITE(4,5040) K 5040 FORMAT('1',20X,'K =',I5,5X,'I=HORIZONTAL,J=DOWN',/) C DO 7015 J=1,JMAX LS= (K-1)*JIM + (J-1)*IM + 1 LE=LS-1+IMAX L=1 DO 7014 I=LS,LE N=ICL(I) IF(N.EQ.-5) THEN NC(L)='-5' ELSE IF(N.EQ.-4) THEN NC(L)='**' ELSE IF(N.EQ.-3) THEN NC(L)='**' ELSE IF(N.EQ.-2) THEN NC(L)='-2' ELSE IF(N.EQ.-1) THEN NC(L)='-1' ELSE IF(N.EQ.0) THEN NC(L)=' 0' ELSE IF(N.GT.0) THEN NC(L)='**' END IF L=L+1 7014 CONTINUE WRITE(4,5030)(NC(L),L=1,IMAX) 5030 FORMAT(30(1X,A2),/) 7015 CONTINUE C 7020 CONTINUE REWIND 4 RETURN END C C C C ================================= SUBROUTINE VERTEX(JA,JB,JC,JVERT) C ================================= C C---- takes the tetrahedron formed by the atom ao and the three C atoms in box(i,j) numbered j=ja,jb,jc. computes the C intersection of the 3 planes between ja,jb,jc,and ao. C this vertex is stored in the array polyhn(i,jvert). C i=1-3=xyz, i=4-6=ja,jb,jc. jvert=serial number of calculated C vertex. C PARAMETER (MAXATM=40000) PARAMETER (MAXSIZ=80000) LOGICAL VRTATM INTEGER SEQNUM,AOSEQ,AONUM CHARACTER*4 ATM,RES3,ATNAMA,RESTYP COMMON /CHRVOL/ ATM(MAXATM),RES3(MAXATM),ATNAMA(200), + RESTYP(30) COMMON/COORDS/ + KEY(MAXATM),NUMCHN(MAXATM),NUMFIL(MAXATM), 1 SEQNUM(MAXATM),X(MAXATM), 2 Y(MAXATM),Z(MAXATM),RVDW(MAXATM),RCOV(MAXATM), 3 LSTATM,LSTRES,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX COMMON/LATICE/ ICL(MAXSIZ), 1 JSERN(MAXATM*2),IMAX,JMAX,KMAX,CE,CE2 COMMON/INPOUP/ ARANGE,VOLUME,METHOD,RDELTA COMMON/SOLVNT/ MAXSOL,NSOL,NMIN2,NTOT,NPAT COMMON/ATDATA/ INDX1(30),INDX2(30),IBND1(200), + KATNAMA(200) COMMON/ORIGN/ AO(3),AOCOV,AOVDW,AOVDW2,AOSEQ,AONUM, + BNDCOD,SEQDIF COMMON/ATMBOX/ NACUBE,NASPHR,IBOX(2,100),VRTATM(100), + VERVOL(100), + BOX(5,100),TEMP(4,100),IPOLYH(4,100),POLYHN(3,100) COMMON/FMRTST/NCALLS(25),JLEVEL,IMSG,IERROR(2,5),ICHK,JCHK DIMENSION PLANES(4,3),P(3) SAVE C C---- Initialize error variables C 80 JLEVEL=JLEVEL+1 DO 90 J=JLEVEL,5 DO 90 I=1,2 90 IERROR(I,J)=0 IERROR(1,JLEVEL)=14 C C---- load array planes(i,j) with data from the array box C concerning ao and the atoms identified by ja,jb,jc. C DO 100 I=1,4 PLANES(I,1)=BOX(I+1,JA) PLANES(I,2)=BOX(I+1,JB) 100 PLANES(I,3)=BOX(I+1,JC) C C---- calculate the intersection of planes C 1450 IERROR(2,JLEVEL)=1450 CALL POINT(PLANES,P) IF(ICHK) 1455,1510,1510 1455 WRITE(3,1) 1 FORMAT(' ERROR IN POINT.VERTEX INDETERMINATE.', + ' RETURN TO ATVOL',/) GO TO 1610 C 1510 DO 1520 I=1,3 1520 POLYHN(I,JVERT)=P(I) IPOLYH(1,JVERT)=JA IPOLYH(2,JVERT)=JB IPOLYH(3,JVERT)=JC 1610 JLEVEL=JLEVEL-1 5000 RETURN END C C C C =================== FUNCTION VOLTET (A) C =================== C C---- returns the volume of the tetrahedron defined by the vertices C given in the array a(3,4). a(i=1,3,j) are the x,y,z C coordinates of the point j. C DIMENSION A(3,4) C CALL DET3C(A,B1,B2,B3,B4) VOLTET= ABS((-B1+B2-B3+B4)/6.) END