hbdata.f

最经典的分子对结软件

      PROGRAM GRIDCONT

****************************************************************
* program to control the grid calculations for all the required
* probes - calculates grid for each probr in turn, then renames
* files as required
****************************************************************

      IMPLICIT NONE

      INTEGER POINTER(1000000) ! used to relate ranked energy to coords
      INTEGER MAXPROBE ! no. of probe atoms keyed for use
      PARAMETER (MAXPROBE = 14)

      INTEGER NUMPLANE ! no. of planes / ang of grid
      INTEGER NUMPROBE ! no. of probes to be processed during current run
      INTEGER * 2 COMPARE ! sorting function
      INTEGER COUNT0,COUNT,COUNT1,COUNT2,K,END,END1,CHECK,NUMEN ! counters
      INTEGER NUMMIN ! used to count no. of minima found
      INTEGER ICRUD ! redundant integer data
      INTEGER  NUMX, NUMY, NUMZ ! no. of x/y/z planes in map
      INTEGER X,Y,Z ! energy array counters
      INTEGER XTEST,YTEST,ZTEST ! used to check energies of gridpoints
*                                 around current test point
      INTEGER NAMELENGTH(MAXPROBE) ! no. of character in probe name - used to
*                            ! to create grid control deck input name

      REAL RCRUD ! redundant real data
      REAL GRIDINC ! grid increment used in map
      REAL MINX, MINY, MINZ ! min x/y/z coords
      REAL EMAX ! max. acceptable energy
      REAL ENERGY(100,100,100),ENERGY1D(4,1000000) ! energy data from map
      REAL XCOORD,YCOORD,ZCOORD ! coord data for energy minima
      REAL MINIMA(4,10000) ! minima coords and energies
      REAL MINRANK(2,100) ! minima ranking variable
      REAL EXCLUDE ! min distance allowed between minima
      REAL XLOW,YLOW,ZLOW,XHIGH,YHIGH,ZHIGH ! grid bounds
      REAL DIST ! distance between energy coords

      CHARACTER * 72 CCRUD ! redundant character data
      CHARACTER * 80 MAPNAME
      CHARACTER * 3 TYPE ! sybyl atom type of probe
      CHARACTER * 3 PROBE(MAXPROBE) ! probe names for grid
      CHARACTER * 80 PROTNAME ! protein name 
      CHARACTER* 80 INDATA,OUTDATA,INCNT!grid calculation control file in/output
	  CHARACTER * 80 DIR,GRINOUT ! dir / filename specs
	  CHARACTER * 80 DATA ! GRIN output listin strings

      EXTERNAL COMPARE

      COMMON /EDATA/ ENERGY1D

      DATA PROBE /"n3+","o::","n1","n2","o","oh","oc2","on","n:#",
     &            "n:=","o=","po4","s1","o1"/

*      DATA PROBE /"n2:","o"/

      DATA NAMELENGTH /3,3,2,2,1,2,3,2,3,3,2,3,2,2/

***************************************************
* first open and read the name of the protein
* and the required grid boundaries
***************************************************

      OPEN (1,STATUS='OLD',FILE="HBIN")

      READ(1,'(A80)') PROTNAME
      READ(1,'(A80)') DIR
      READ(1,'(F8.3)') XLOW
      READ(1,'(F8.3)') YLOW
      READ(1,'(F8.3)') ZLOW
      READ(1,'(F8.3)') XHIGH
      READ(1,'(F8.3)') YHIGH
      READ(1,'(F8.3)') ZHIGH
      READ(1,'(I1)') NUMPLANE
      READ(1,'(F6.3)') EMAX
      READ(1,'(F5.3)') EXCLUDE
      READ(1,'(I2)') NUMPROBE 
     
      CLOSE (1) 

*************************************
* extract protein name from string
*************************************

      END = 1
1000  IF (PROTNAME(END:END).EQ.' ') THEN
         GOTO 2000
      ELSE
         END = END + 1
         GOTO 1000
      END IF

*************************************
* extract grid directory from string
*************************************


2000      END1 = 1
3000  IF (DIR(END1:END1).EQ.' ') THEN
         GOTO 4000
      ELSE
         END1 = END1 + 1
         GOTO 3000
      END IF



**********************************
* set up GRIN file fro pdb prep
**********************************

4000  OPEN(1,STATUS='UNKNOWN',FILE='GRININ')


	  WRITE(1,'(A1)') '1'
	  WRITE(1,'(A1)') '1'
	  WRITE(1,'(A1)') '1'
	  WRITE(1,'(A80)') PROTNAME(1:END-1) // '.pdb'
	  WRITE(1,'(A1)') '2' 
	  WRITE(1,'(A80)') PROTNAME(1:END-1) // '_grin.lis' 
	  WRITE(1,'(A1)') '3'
	  WRITE(1,'(A80)') PROTNAME(1:END-1) // '_grin.dat'   
	  WRITE(1,'(A1)') '4'
	  WRITE(1,'(A80)') DIR(1:END1-1) // '/grin' 
	  WRITE(1,'(A1)') '6'
	  WRITE(1,'(A80)') DIR(1:END1-1) // '/grub.dat' 
	  WRITE(1,'(A1)') '9'
	  WRITE(1,'(A1)') '4'
	  WRITE(1,'(A2)') 'no'
	  WRITE(1,'(A2)') 'no'
	  WRITE(1,'(A7)') 'certain'
	  WRITE(1,'(A1)') ''
	  WRITE(1,'(A1)') '7'
	  WRITE(1,'(A1)') '7'
	  WRITE(1,'(A1)') '7'
	  WRITE(1,'(A1)') '7'

	  CLOSE(1)


*************
* run GRIN
*************

      DATA = DIR(1:END1-1) // '/great < GRININ'


      CALL SYSTEM (DATA)


*******************************
* scan out listing for errors
*******************************

       GRINOUT=PROTNAME(1:END-1) // '_grin.lis'

	   OPEN (1,STATUS='OLD',FILE=GRINOUT)

5000   READ(1,'(A80)',end=6000) DATA
	   IF (DATA(8:55) .EQ. 
     & 'BUT IT IS IMPORTANT TO CHECK YOUR PDB INPUT FILE') THEN
	        WRITE(*,*) ' probable error in pdb file - check ',
     &      GRINOUT(1:END+10), 'for details'
   	     	STOP
	   END IF

       GOTO 5000

6000   CLOSE(1)

**************************************************************
* for each probe type create the required input deck
**************************************************************
     
	  DO  COUNT = 1,NUMPROBE
 
        OPEN(1,STATUS='UNKNOWN',FILE='GRIDIN') 
   

**********************************************************
* first write template info which tells grid where where
* its executables are
**********************************************************

 
        WRITE(1,'(A1)') '8'
        WRITE(1,'(A1)') '2'
        WRITE(1,'(A1)') '1'
        WRITE(1,'()')
         

****************************************
* write out required output filenames
****************************************

	    WRITE(1,'(A1)') '1'
	    WRITE(1,'(A80)') PROTNAME(1:END-1) // '_grin.dat'   
        WRITE(1,'(A1)') '2'
        INDATA = PROTNAME(1:(END-1)) // '_' // 
     &            PROBE(COUNT)(1:NAMELENGTH(COUNT)) // '.lis'
        WRITE(1,'(A80)') INDATA

        WRITE(1,'(A1)') '3'
        INDATA = PROTNAME(1:(END-1)) // '_' // 
     &            PROBE(COUNT)(1:NAMELENGTH(COUNT)) // '.dat'
        WRITE(1,'(A80)') INDATA

        WRITE(1,'(A1)') '4'
        WRITE(1,'(A80)') '/bert/grid11.01/grid'

        WRITE(1,'(A1)') '6'
        INDATA = PROTNAME(1:(END-1)) // '_' // 
     &            PROBE(COUNT)(1:NAMELENGTH(COUNT)) // '.in'
        WRITE(1,'(A80)') INDATA
       

*************************************
* set required grid boundaries
*************************************
 
        WRITE(1,'(A1)') '8'
        WRITE(1,'(A1)') '4'
        WRITE(1,'(A1)') '2'
        WRITE(1,'(3(F8.3,1X))') XHIGH,YHIGH,ZHIGH
        WRITE(1,'(A1)') 'Y'
        WRITE(1,'(3(F8.3,1X))') XLOW,YLOW,ZLOW
        WRITE(1,'(A1)') 'Y'


***************************
* set up grid density
***************************

        WRITE(1,'(A1)') '8'
        WRITE(1,'(A1)') '8'
        WRITE(1,'(A1)') '1'
        WRITE(1,'(A1)') '6'
        WRITE(1,'(I1)') NUMPLANE
         

****************************
* write out required probe
****************************

        WRITE(1,'(A1)') '8'
        WRITE(1,'(A1)') '8'
        WRITE(1,'(A1)') '5'
        WRITE(1,'(A1)') '3'
        WRITE(1,'(A3)') PROBE(COUNT)
        WRITE(1,'(A1)') 'Y'


***************************************************
* write execution, comment and terminate commands
***************************************************

        WRITE(1,'(A1)') '9'

        WRITE(1,'(A1)') '4'
        INDATA = PROTNAME(1:(END-1)) // ' with ' // PROBE(COUNT) //
     &            ' calculation'
        WRITE(1,'(A80)') INDATA

        WRITE(1,'(A1)') 'Y'
        WRITE(1,'()')
        WRITE(1,'(A1)') '7'
        WRITE(1,'(A1)') '7'
        WRITE(1,'(A1)') '7'
        WRITE(1,'(A1)') '7'
      
        CLOSE (1)


*************************************************
* spawn command to system so as to run grid
*************************************************

		DATA = DIR(1:END1-1) // '/great < GRIDIN'

        CALL SYSTEM (DATA)

        INCNT = PROTNAME(1:(END-1)) // '_' //
     &            PROBE(COUNT)(1:NAMELENGTH(COUNT)) // '.cnt'

        INDATA = PROTNAME(1:(END-1)) // '_' //
     &            PROBE(COUNT)(1:NAMELENGTH(COUNT)) // '.dat'

        OUTDATA = PROTNAME(1:(END-1)) // '_' //
     &          PROBE(COUNT)(1:NAMELENGTH(COUNT)) // '.mol2'

        OPEN (10,STATUS='OLD',FORM='UNFORMATTED',FILE=INDATA)


****************************************************************
* read header data - most of it is redundant - extract grid
* data only
****************************************************************

        READ(10) CCRUD,RCRUD,RCRUD,ICRUD,ICRUD,ICRUD,ICRUD,ICRUD,
     &           NUMX,NUMY,NUMZ,GRIDINC,MINX,MINY,MINZ


**********************************************************
* use header data to extract energy data from map
**********************************************************

        DO Z = 1,NUMZ
          READ(10)
          READ(10) ((ENERGY(X,Y,Z), X=1,NUMX), Y=1,NUMY)
        END DO

        NUMEN = 0
   
        DO Z=1,NUMZ
           DO Y = 1,NUMY
              DO X = 1,NUMX

                NUMEN = NUMEN + 1
                ENERGY1D(1,NUMEN) = MINX + (GRIDINC*X)
                ENERGY1D(2,NUMEN) = MINY + (GRIDINC*Y)
                ENERGY1D(3,NUMEN) = MINZ + (GRIDINC*Z)
                ENERGY1D(4,NUMEN) = ENERGY(X,Y,Z)

              END DO
           END DO
        END DO


        DO X = 1, NUMX*NUMY*NUMZ
           POINTER(X) = X
        END DO

*****************
* sort energies 
*****************

      CALL QSORT(POINTER,100000,4,COMPARE)



***************************************************************
* go down energy list - keep each new lowest energy .gt.
* given mindist from all other minimum energies found so far 
***************************************************************

       NUMMIN = 0
       NUMEN = 0

       DO X = 1,NUMX*NUMY*NUMZ

         NUMEN = NUMEN + 1


*******************************************
* skip over all high energy points
*******************************************

         IF (ENERGY1D(4,POINTER(NUMEN)) .GT. EMAX) GOTO 7000

         DO CHECK = 1, NUMMIN

           DIST=SQRT(((MINIMA(1,CHECK)-ENERGY1D(1,POINTER(NUMEN)))**2)+
     &               ((MINIMA(2,CHECK)-ENERGY1D(2,POINTER(NUMEN)))**2)+
     &               ((MINIMA(3,CHECK)-ENERGY1D(3,POINTER(NUMEN)))**2))                
           IF (DIST .LT. EXCLUDE) GOTO 7000

         END DO     

              NUMMIN = NUMMIN + 1

              MINIMA(1,NUMMIN) = ENERGY1D(1,POINTER(NUMEN)) 
              MINIMA(2,NUMMIN) = ENERGY1D(2,POINTER(NUMEN)) 
              MINIMA(3,NUMMIN) = ENERGY1D(3,POINTER(NUMEN)) 
              MINIMA(4,NUMMIN) = ENERGY1D(4,POINTER(NUMEN))

7000   END DO
   
        
********************************************
* write out minima results into mol2 file
********************************************

      CALL MAKEMOL2(NUMMIN,MINIMA,OUTDATA)


        CLOSE(10)

*****************************************************************
* create input file and use it to create sybyl contour file
*****************************************************************

        OPEN(20,STATUS='UNKNOWN',FILE='GSYBIN')
        WRITE(20,'(A80)') INDATA
        WRITE(20,'(A80)') INCNT
        CLOSE(20)
 
		DATA = DIR(1:END1-1) // '/gsyb < GSYBIN'

        CALL SYSTEM (DATA)

      END DO


      END

*===============================================================

      INTEGER * 2 FUNCTION COMPARE (ONE,TWO)

**************************************
* function for sorting routine QSORT
**************************************


      IMPLICIT NONE

      INTEGER ONE,TWO
      REAL ENERGY1D(4,1000000)

      COMMON /EDATA/ ENERGY1D

      IF (ENERGY1D(4,ONE).GT.ENERGY1D(4,TWO)) THEN
         COMPARE = 1
      ELSE IF (ENERGY1D(4,ONE).LT.ENERGY1D(4,TWO)) THEN
         COMPARE = -1
      ELSE
         COMPARE = 0
      END IF

      RETURN

      END


*==================================================================

      SUBROUTINE MAKEMOL2(NUMMIN,MINIMA,OUTMOL2)

****************************************************
* routine for the creation of mol2 surface files
****************************************************

      IMPLICIT NONE

      INTEGER  COUNT,K ! counters
      INTEGER  NUMMIN ! no. of site points found so far

      REAL MINIMA(4,10000) ! sitepoint data
      REAL GRIDINC
      REAL MAPLIMITS(6) ! min x,y,x / max x,y,z values for grid

      CHARACTER * 80 OUTMOL2 ! name of output mol2 file
      CHARACTER * 3 TYPE ! used to control atom type assigned to sitepoint

      OPEN(20,STATUS='UNKNOWN',FILE=OUTMOL2)

      WRITE(20,'(A17)') '@<TRIPOS>MOLECULE'


***************************************************
* have to write grid data in name field for lf
***************************************************

      WRITE(20,'(A80)') OUTMOL2
      WRITE(20,'(5I4)') NUMMIN,0,0,0,0
      WRITE(20,'(A5)') 'SMALL'
      WRITE(20,'(A12)') 'USER_CHARGES'
      WRITE(20,'()')
      WRITE(20,'()')
      WRITE(20,'(A13)') '@<TRIPOS>ATOM'


      DO COUNT = 1, NUMMIN

        WRITE(20,'(I4,1X,A1,I4.4,1X,3F9.4,1X,A3,1X,I1,1X,A3,1X,F6.2)')
     &  COUNT,"S",COUNT,MINIMA(1,COUNT),MINIMA(2,COUNT),
     &  MINIMA(3,COUNT),"du ",0,"<0>",MINIMA(4,COUNT)

      END DO

      CLOSE (20)

      RETURN

      END