C C $Id: clustsub.F,v 1.11 1998/07/16 16:39:31 jjv5 Exp $ C C------------------------------------------------------------------------ SUBROUTINE MAKESUB(ipoint, natot, nbig, ibig, IERROR) C C ROUTINE TO DO CLUSTER SUBSETTING. C C March 31 1998, Arjan van der Vaart: added multiple core feature. C the core will consist of ncore(1) residues selected from the C group of residues icorel(1)..icorel(icorel1(2)-1), ncore(2) residues C selected form the group of residues C icorel(icorel1(2))..icorel(icorel1(3)-1) etc. Buffers will always C be build from the entire set (1..nres) of residues. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" C C C LOCAL: C DIMENSION IRCORE(MAXRES),IR1(MAXRES), . ILOOK(MAXRES),ICATMS(MAXATM),IBATMS(2,MAXATM), & ibuflist(maxres), itmp(maxres), rijlist(MAXATM), iatmlist(MAXATM) LOGICAL INCORE(MAXRES),INSUB(MAXRES),INGSUB(MAXRES),HARSH, & addineighb C CC-----------------------------------------------------------------------CC #ifdef MEMORY_OVERLAP C now the following arrays share memory: C array1 common | array2 common | shared memory array C ----------------------------------------------------------------------- C rdist(mxpres) /clst/ | ff(msorb2) /work/ | ff(maxovlp1) C #define RDIST_ ff #else C no sharing of memory of arrays #define RDIST_ rdist #endif CC---------------------------------------------------------------------CC addineighb = .false. if ((.not.noovrlp).and.(nneighbor.gt.0)) then C-RDC write(iout,'(/" WARNING: IGNORED NEIGHBOR PARAMETERS SINCE ", C-RDC & "NO-OVERLAP KEYWORD ", C-RDC & /" IS NOT DEFINED IN divcon.in")') nneighbor = 0 endif if (nneighbor.gt.0) then if (neighbor(nneighbor).gt.ncores) then C-RDC write(iout,'(/" ERROR: NEIGHBOR ",I5," IS OUT OF RANGE", C-RDC & /" NEIGHBORS SHOULD BE BETWEEN 1 AND ",I5)') C-RDC & neighbor(nneighbor), ncores ierror = 1 return endif do i=1,nneighbor if (ncore(neighbor(i)).ne.1) then C-RDC write(iout,'(/" ERROR: NEIGHBOR GROUP ",I5," DOES NOT ", C-RDC & /" HAVE ONE RESIDUE IN THE CORE")') C-RDC & neighbor(i) ierror = 1 return endif enddo neighbor(nneighbor+1) = 0 else neighbor(1) = 0 endif ineighb = 1 IF(PRTSUB) WRITE(IOUT,'(" SUBSYSTEM DEFINITIONS:")') C C RCUT = SQUARE OF ACTUAL CUTOFF IN ANGSTROMS. RCUT = 49.0D0 c RBUFF1 = DCBUFF1**2 c RBUFF2 = (DCBUFF1 + DCBUFF2)**2 NLARGE = 0 ISUBMX = 0 ncsub = 0 ncatot = 0 isubend(0) = 0 C make cores for all different values of ncore do 5000 icr=1,ncores if (addineighb) then ineighb = ineighb + 1 addineighb = .false. endif C . i0 is the offset for the residue list and nresn is the C . number of residues with ncore=ncore(icr) i0 = icorel1(icr)-1 nresn = icorel1(icr+1)-icorel1(icr) IF(NCORE(ICR).NE.1.AND.NCORE(ICR).NE.NRESN)THEN C C . CONSTRUCT A RESIDUE-BASED PAIRLIST AND STORE THE MININUM C . INTER-RESIDUE INTERATOMIC DISTANCE IN EACH CASE. STORE C . PAIRS THAT ARE WITHIN 5.0 ANGSTROMS. IF THERE ARE ISOLATED C . RESIDUES THAT ARE LOCATED A GREAT DISTANCE FROM ALL OTHER C . RESIDUES, THEN IT MAY BE NECESSARY TO INCREASE THIS CUTOFF C . IN ORDER TO GET SUCCESSFUL SUBSETTING. C ISTORE = 1 DO 80 IR=1,NRESN ires = icorel(ir+i0) IR1(IR) = ISTORE IA1 = IRPNT(IRES) IA2 = IRPNT(IRES+1)-1 DO 70 JR=1,NRESN jres = icorel(jr+i0) IF(IRES.EQ.JRES) GO TO 70 JA1 = IRPNT(JRES) JA2 = IRPNT(JRES+1)-1 RMIN = 1.0D10 DO 60 IATM=IA1,IA2 IF(NATORB(IATNUM(IATM)).EQ.0) GO TO 60 DO 50 JATM=JA1,JA2 IF(NATORB(IATNUM(JATM)).EQ.0) GO TO 50 xj = xyz(1,jatm) yj = xyz(2,jatm) zj = xyz(3,jatm) RSQR = (XYZ(1,IATM)-XJ)**2 . + (XYZ(2,IATM)-YJ)**2 . + (XYZ(3,IATM)-ZJ)**2 RMIN = MIN(RMIN,RSQR) 50 CONTINUE 60 CONTINUE IF(RMIN.LT.RCUT)THEN IRPAIR(ISTORE) = JR RDIST_(ISTORE) = RMIN ISTORE = ISTORE + 1 IF(ISTORE.GT.MXPRES)THEN IERROR = 1 WRITE(IOUT,'(/" MAXIMUM STORAGE FOR", . " RESIDUE-BASED", . " PAIRLIST EXCEEDED",/" -- INCREASE", . " MXPRES PARAMETER IN divcon.dim AND", . " RECOMPILE")') RETURN ENDIF ENDIF 70 CONTINUE C C . SORT PAIRLIST ENTRIES FOR IRES ON THE BASIS OF INCREASING C . DISTANCE. C I1 = IR1(IR) NVAL = ISTORE - I1 CALL RSORT(NVAL,RDIST_(I1),IRPAIR(I1)) 80 CONTINUE IR1(NRESN+1) = ISTORE ENDIF C C . BUILD EACH SUBSYSTEM BY GENERATING A CORE COLLECTION OF RESIDUES, C . THEN EXTENDING OUT USING THE BUFFER THICKNESS. C C . PBCs WILL BE USED (IF REQUESTED) TO GENERATE BUFFER REGIONS. C C NUSED = 0 ISEED1 = 1 DO 90 IRES=1,NRESN INCORE(IRES) = .FALSE. 90 CONTINUE C C . NUSED = COUNT OF RESIDUES THAT APPEAR IN ANY CORE. C C . ISEED1 = STARTING POINT FOR SEARCH IN RESIDUE LIST TO FIND C . SEED RESIDUE TO START NEXT CORE. C C . INCORE = FLAG FOR RESIDUES THAT APPEAR IN ANY CORE. C C . here all arrays and iseed is used locally i.e. filled in from C . 1 till nresn. This means that e.g. incore(i) means that residue C . icorel(i+i0) is in the core. C . ingsub(i) is the only global array i.e. the indices go from 1 C . till nres. This means that ingsub(i) means that residue i is C . included in the current subsystem. C DO 500 IS=1,NRESN isub = icorel(is+i0) nsub = nsub + 1 ncsub = ncsub + 1 HARSH = .TRUE. C C . INITIALLY, A HARSH CRITERION WILL BE USED FOR ADDITION OF C . RESIDUES TO CORE ONCE THE SEED RESIDUE HAS BEEN SELECTED. C . THE HARSH CRITERION IS THAT A RESIDUE WILL NOT BE ADDED C . TO THE CORE OF ISUB IF IT APPEARS IN THE CORE OF ANY OTHER C . SUBSYSTEM. IF RESIDUES BECOME "ISOLATED" OR IF NCORE DOES C . NOT DIVIDE EVENLY INTO NRES, THEN IT MAY BECOME NECESSARY TO C . RELAX CRITERION AND ALLOW PREVIOUSLY USED RESIDUES TO BE ADDED. C 100 NADDED = 0 ISEED0 = ISEED1 C C . NADDED = THE NUMBER OF RESIDUES IN THE CORE OF ISUB. C C . ISEED0 WILL BE USED TO RESET ISEED1 IF THE SEARCH HAS TO C . BE REPEATED BECAUSE THE HARSH CRITERION WAS TOO HARSH. C DO 110 IRES=1,NRESN INSUB(IRES) = .FALSE. 110 CONTINUE do ires=1,nres ingsub(ires) = .false. enddo C C . INSUB FLAGS RESIDUES THAT HAVE BEEN ADDED TO THE CORE OF ISUB. C IF(NCORE(ICR).EQ.1)THEN C C . ONLY ONE RESIDUE PER CORE -- ISUB WILL JUST CORRESPOND TO C . THE CORE RESIDUE. C NUSED = NUSED + 1 NADDED = NADDED + 1 IRCORE(1) = IS INSUB(IS) = .TRUE. ingsub(isub) = .true. ELSEIF(NCORE(ICR).EQ.NRESN)THEN C C . ONLY ONE SUBSYSTEM -- ALL RESIDUES. C NUSED = NRESN NADDED = NRESN DO 120 IRES=1,NRESN IRCORE(IRES) = IRES ingsub(icorel(ires+i0)) = .true. 120 CONTINUE ELSE C C . SEARCH RESIDUES THAT DO NOT ALREADY APPEAR IN ANY CORE LIST TO C . FIND A SEED AS A STARTING POINT FOR THE NEXT CORE. C DO 140 ISEED=ISEED1,NRESN IF(.NOT.INCORE(ISEED))THEN INCORE(ISEED) = .TRUE. INSUB(ISEED) = .TRUE. ingsub(icorel(iseed+i0)) = .true. ILOOK(ISEED) = IR1(ISEED) C C . ILOOK(IRES) = STARTING POINT FOR SEARCH IN IRES PAIRLIST C . TO LOCATE THE NEXT "CLOSE" RESIDUE. C NUSED = NUSED + 1 NADDED = 1 IRCORE(1) = ISEED C C . IRCORE(I) = ITH CORE RESIDUE IN ISUB. C GO TO 150 ENDIF 140 CONTINUE 150 ISEED1 = ISEED + 1 C C . NOW MAKE NCORE-1 PASSES THROUGH THE PAIRLIST OF THE CORE C . RESIDUES, BRINGING IN THE CLOSEST RESIDUE JRMIN AT EACH PASS. C DO 200 IPASS=2,NCORE(ICR) RMIN = 1.0D10 JRMIN = 0 DO 180 ICORE=1,IPASS-1 IR = IRCORE(ICORE) ires = icorel(ir+i0) C C . SKIP THIS CORE RESIDUE IF WE'VE REACHED THE END OF C . ITS PAIRLIST. C 160 J1 = ILOOK(IR) J2 = IR1(IR+1) IF(J1.EQ.J2) GO TO 180 ILOOK(IR) = J1+1 C JRES = IRPAIR(J1) C C . IF JRES HAS ALREADY BEEN USED (WHICHEVER CRITERION), C . THEN GO BACK AND CHECK THE NEXT ENTRY IN THE IRES C . PAIRLIST. C IF(HARSH)THEN IF(INCORE(JRES)) GO TO 160 ELSE IF(INSUB(JRES)) GO TO 160 ENDIF C RIJ = RDIST_(J1) IF(RIJ.LT.RMIN)THEN RMIN = RIJ JRMIN = JRES ENDIF 180 CONTINUE IF(JRMIN.EQ.0.AND..NOT.HARSH)THEN C C . RETURN WITH AN ERROR IF THE RELAXED CRITERION IS C . IN EFFECT AND WE STILL CAN'T FIND A CLOSE RESIDUE C . IN THE CORE PAIRLIST. IT IS PROBABLY A RESULT OF C . THE CURRENT CORES BEING FAR AWAY FROM EVERYTHING ELSE. C . IT SHOULD BE POSSIBLE TO RESOLVE THIS PROBLEM BY C . INCREASING RCUT. C IERROR = 1 WRITE(IOUT,'(/" CANNOT BUILD CORE FOR SUBSYSTEM ", . I4,/" -- TRY INCREASING CUTOFF VALUE IN", . " CLUSTSUB")') ISUB RETURN ELSEIF(JRMIN.EQ.0)THEN C C . COULDN'T FIND A CLOSE RESIDUE USING HARSH CRITERION. C . SWITCH TO RELAXED CRITERION AND TRY AGAIN. C HARSH = .FALSE. DO 190 I=1,IPASS-1 IRES = IRCORE(I) INCORE(IRES) = .FALSE. 190 CONTINUE NUSED = NUSED - NADDED ISEED1 = ISEED0 GO TO 100 ENDIF C C . ADD CLOSEST RESIDUE TO CORE. C IRCORE(IPASS) = JRMIN ILOOK(JRMIN) = IR1(JRMIN) INSUB(JRMIN) = .TRUE. ingsub(icorel(jrmin+i0)) = .true. IF(HARSH.AND.INCORE(JRMIN)) GO TO 200 IF(.NOT.INCORE(JRMIN)) NUSED = NUSED + 1 NADDED = NADDED + 1 INCORE(JRMIN) = .TRUE. 200 CONTINUE ENDIF C C . CORE IS COMPLETE. STORE CORE ATOMS IN ICATMS. C NACORE = 0 DO 220 IADD=1,NADDED IR = IRCORE(IADD) ires = icorel(ir+i0) IA1 = IRPNT(IRES) IA2 = IRPNT(IRES+1)-1 DO 210 IATM=IA1,IA2 IF(NATORB(IATNUM(IATM)).EQ.0) GO TO 210 NACORE = NACORE + 1 ICATMS(NACORE) = IATM 210 CONTINUE 220 CONTINUE C C . CREATE TWO SETS OF BUFFER ATOMS FOR THIS SUBSYSTEM. THE C . FIRST SET WILL CONTAIN ALL ATOMS THAT ARE WITHIN SQRT(RBUFF1) C . OF THE CORE. THE SECOND SET WILL CONTAIN ATOMS THAT ARE C . BEYOND SQRT(RBUFF1) BUT WITHIN SQRT(RBUFF2) OF CORE. C C . if noovrlp is .false. then C . note that atoms for the buffers are extracted from the global C . set of residues, i.e. from residues 1..nres (and not only C . from residues icorel(1+i0)..icorel(nresn+i0)) C C . if noovrlp is .true. then C . buffers are created from "own" local list C NBUFF1 = 0 NBUFF2 = 0 ncandidate =0 IF(NCORE(ICR).EQ.NRES) GO TO 400 if (.not.noovrlp) then DO 300 IRES=1,NRES IF(INGSUB(IRES)) GO TO 300 IA1 = IRPNT(IRES) IA2 = IRPNT(IRES+1)-1 DO 280 IATM=IA1,IA2 C C . IATM IS THE CANDIDATE ATOM FOR ADDITION TO BUFFER LISTS. C IF(NATORB(IATNUM(IATM)).EQ.0) GO TO 280 XI = XYZ(1,IATM) YI = XYZ(2,IATM) ZI = XYZ(3,IATM) C C . GET MINIMUM DISTANCE BETWEEN IATM AND CORE ATOMS. C RMIN = 1.0D10 DO 250 J=1,NACORE JATM = ICATMS(J) IF(PBC)THEN CALL PBCXYZ(IATM,JATM,XJ,YJ,ZJ) ELSE XJ = XYZ(1,JATM) YJ = XYZ(2,JATM) ZJ = XYZ(3,JATM) ENDIF RIJ = (XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2 RMIN = MIN(RMIN,RIJ) 250 CONTINUE ncandidate = ncandidate + 1 iatmlist(ncandidate)=iatm rijlist(ncandidate)=rmin cc IF(RMIN.LE.RBUFF1)THEN cc NBUFF1 = NBUFF1 + 1 cc IBATMS(1,NBUFF1) = IATM cc ELSEIF(RMIN.LE.RBUFF2)THEN cc NBUFF2 = NBUFF2 + 1 cc IBATMS(2,NBUFF2) = IATM cc ENDIF 280 CONTINUE 300 CONTINUE do i=1,natoms if (i.gt.ncandidate) rijlist(i)=1.0d10 enddo call rsort(natoms,rijlist,iatmlist) do i=1,natbuff1 IBATMS(1,i) = iatmlist(i) enddo do i=1,natbuff2 IBATMS(2,i) = iatmlist(i+natbuff1) enddo else if (neighbor(ineighb).eq.icr) then ir = ircore(1) ires = icorel(i0+ir) ileft = ires - neighborn(ineighb) ileft = max(1,ileft) ntmp = 0 do j=ileft,ires-1 ntmp = ntmp + 1 itmp(ntmp) = j enddo iright = ires + neighborn(ineighb) iright = min(nres,iright) do j=ires+1,iright ntmp = ntmp + 1 itmp(ntmp) = j enddo j0 = 0 if (ntmp.gt.0) then k = 1 do j=1,nresn 301 if (icorel(j+i0).gt.itmp(k)) then k = k + 1 if (k.gt.ntmp) goto 302 goto 301 elseif (icorel(j+i0).eq.itmp(k)) then j0 = j0 + 1 ibuflist(j0) = itmp(k) k = k + 1 if (k.gt.ntmp) goto 302 endif enddo endif 302 addineighb = .true. else j0 = nresn do j=1,nresn ibuflist(j) = icorel(j+i0) enddo endif C the following commented code is without neighbors c do 390 ir=1,nresn C ires = icorel(ir+i0) ncandidate =0 do 390 ir=1,j0 ires = ibuflist(ir) IF(INGSUB(IRES)) GO TO 390 IA1 = IRPNT(IRES) IA2 = IRPNT(IRES+1)-1 DO 380 IATM=IA1,IA2 C C . IATM IS THE CANDIDATE ATOM FOR ADDITION TO BUFFER LISTS. C IF(NATORB(IATNUM(IATM)).EQ.0) GO TO 380 XI = XYZ(1,IATM) YI = XYZ(2,IATM) ZI = XYZ(3,IATM) C C . GET MINIMUM DISTANCE BETWEEN IATM AND CORE ATOMS. C RMIN = 1.0D10 DO 350 J=1,NACORE JATM = ICATMS(J) IF(PBC)THEN CALL PBCXYZ(IATM,JATM,XJ,YJ,ZJ) ELSE XJ = XYZ(1,JATM) YJ = XYZ(2,JATM) ZJ = XYZ(3,JATM) ENDIF RIJ = (XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2 RMIN = MIN(RMIN,RIJ) 350 CONTINUE ncandidate = ncandidate + 1 iatmlist(ncandidate)=iatm rijlist(ncandidate)=rmin c IF(RMIN.LE.RBUFF1)THEN c NBUFF1 = NBUFF1 + 1 c IBATMS(1,NBUFF1) = IATM c ELSEIF(RMIN.LE.RBUFF2)THEN c NBUFF2 = NBUFF2 + 1 c IBATMS(2,NBUFF2) = IATM c ENDIF 380 CONTINUE 390 CONTINUE do i=1,natoms if (i.gt.ncandidate) rijlist(i)=1.0d10 enddo call rsort(natoms,rijlist,iatmlist) do i=1,natbuff1 IBATMS(1,i) = iatmlist(i) enddo do i=1,natbuff2 IBATMS(2,i) = iatmlist(i+natbuff1) enddo endif C C-RDC C . WRITE OUT LATEST SUBSYSTEM INFORMATION IF REQUESTED. C 400 IF(PRTSUB)THEN WRITE(IOUT,'(/" SUBSYSTEM ",I4, . /" --------------", . /" CORE ATOMS:")') NSUB WRITE(IOUT,'(16I5)') (ICATMS(I),I=1,NACORE) IF(NatBUFF1.GT.0)THEN WRITE(IOUT,'(" INNER BUFFER LAYER:")') WRITE(IOUT,'(16I5)') (IBATMS(1,I),I=1,NBUFF1) ENDIF IF(NatBUFF2.GT.0)THEN WRITE(IOUT,'(" OUTER BUFFER LAYER:")') WRITE(IOUT,'(16I5)') (IBATMS(2,I),I=1,NBUFF2) ENDIF WRITE(IOUT,*) ENDIF C NAT = NACORE + NatBUFF1 + NatBUFF2 C C . MAKE SURE WE HAVE ENOUGH SPACE TO STORE THE NAT ATOMS IN C . THIS SUBSYSTEM. C IF((IPOINT+NAT).GE.MSLIST)THEN IERROR = 1 WRITE(IOUT,'(/" ATOM LIST STORAGE LIMIT REACHED ", . "AT SUBSYSTEM",i4, . /" -- INCREASE MSLIST PARAMETER IN dicon.dim", . " AND RECOMPILE")') NSUB RETURN ENDIF C C . STORE CORE AND BUFFER ATOMS IN SUBSYSTEM ATOM LIST. ASSIGN C . BUFFER STATUS IABUFF AS WE GO C IATOM1(NSUB) = IPOINT DO 420 I=1,NACORE IATOMS(IPOINT) = ICATMS(I) IABUFF(IPOINT) = 0 IPOINT = IPOINT + 1 420 CONTINUE IF(NatBUFF1.GT.0)THEN DO 440 I=1,NatBUFF1 IATOMS(IPOINT) = IBATMS(1,I) IABUFF(IPOINT) = 1 IPOINT = IPOINT + 1 440 CONTINUE ENDIF IF(NatBUFF2.GT.0)THEN DO 460 I=1,NatBUFF2 IATOMS(IPOINT) = IBATMS(2,I) IABUFF(IPOINT) = 2 IPOINT = IPOINT + 1 460 CONTINUE ENDIF C C . SORT ATOM LIST AND CARRY ALONG BUFFER STATUS. C ISTART = IATOM1(NSUB) CALL BSORT(NAT,IATOMS(ISTART),IABUFF(ISTART)) C NATOT = NATOT + NAT ncatot = ncatot + nat IF(NAT.GT.NLARGE)THEN NLARGE = NAT ISUBMX = NSUB if (nat.gt.nbig) then nbig = nat ibig = nsub endif ENDIF C IF(NUSED.EQ.NRESN)THEN C C . ALL RESIDUES HAVE BEEN USED -- WE ARE DONE SUBSETTING. C goto 4900 endif 500 CONTINUE 4900 isubend(icr) = nsub 5000 continue C OVERFLOW POINTER: C IATOM1(NSUB+1) = IPOINT C NATAVG = NCATOT/NCSUB C-RDC WRITE(IOUT,'(//" SUMMARY OF CLUSTER BASED SUBSETTING:", C-RDC . /" -----------------------------------", C-RDC . //" NUMBER OF SUBSYSTEMS = ",I4, C-RDC . /" AVERAGE SUBSYSTEM SIZE = ",I4, C-RDC . /" MAXIMUM SUBSYSTEM SIZE = ",I4, C-RDC . /" LARGEST SUBSYSTEM = #",I4)') C-RDC . NCSUB,NATAVG,NLARGE,ISUBMX if (screen) then C-RDC write(iscr,'(/"cluster subsetting: nsub = ",i4,", av. size = ", C-RDC & i4," max. size = ",i4)') ncsub,natavg,nlarge endif RETURN END CC--------------------------------------------------------------------CC #undef RDIST_